ランダム効果と
correlation structure

ここでは、ランダム効果の導入が反復測定デザインの解析に有効な理由について、"correlation structure"を用いた直観的な説明を試みる。

Correlation structureとは、ざっくりいうと「データや効果の類似性（の期待値）の構造」のことである。

Rでは、誤差間に特定のcorrelation structureを仮定した解析を行うことができる。

面白いことに、この方法を用いると、ランダム効果を導入することなしに、ランダム効果を導入する方法と全く同じ結論にたどり着く場合がある。このことは、ランダム効果とcorrelation structureとの間の密接な関係を示唆している。

ランダム効果の導入がcorrelation structureにどのような影響を及ぼすのか理解することは、ランダム効果の直観的なイメージを得る助けにもなるはずである。

要点

• Correlation structureは”総当たり戦表（正方行列）”で表現される

• ANOVAや通常の線形回帰モデルでは、各ペアのcorrelation（類似性）はすべてゼロである

• ランダム切片モデルでは、ランダム効果のグループ内のすべてのペアで、correlationが等しくなる

• Correlation structureを直接指定することで、さらに柔軟なモデリングを行うことができる

1. 一般線形モデルのcorrelation structure

1.1. 一般線形モデルの作成

昆虫の幼虫の体重を経時的に3回測定し、以下のようなデータが得られたとする。Weightは体重、Timeは測定時点、Larval_IDは幼虫の個体である。ここで、時間経過に従って体重がどのように増加していくのか知りたいとする。

>dat

     Weight Time Larval_ID

1 0.8625486    1         1

2 1.2923105    2         1

3 1.4263123    3         1

4 1.9185645    1         2

5 2.3913196    2         2

6 2.1835319    3         2

7 1.8324626    1         3

8 1.3362665    2         3

9 1.5063424    3         3

まずは説明のため、個体を考慮しない、単なる一般線形モデルで体重増加を表現してみる。

時間に沿って体重が増加するパターンが表現されている。

係数の値は

• Intercept: 1.54

• Time2: 0.14

• Time3: 0.17

となっている。

また、誤差の標準偏差つまり標準誤差（Residual standard error）は、0.55となっている。

1.2. 一般線形モデルの行列表示

次に、上記のモデルを行列表示することを考える。

イメージがつきやすいよう、一般式を示す前に、まずは具体的な数字を使って行列表示を試みる。

ここで、Yは測定値の縦ベクトル、Xは説明変数の行列（デザイン行列と呼ばれる）、βは係数の縦ベクトル、εは誤差の縦ベクトルである。デザイン行列は、左の列から順に切片、Time2、Time3のダミーコードが並んでいる。またβには上から順に切片、Time2、Time3の計数推定値が並んでいる（上のsummaryと見比べてみよう）。たとえば、一番上の測定値であれば、0.86 = (1×1.54 + 0×0.14 + 0×0.17) + (-0.68) と計算されることがわかる。

では、これを一般式に表してみよう。すると、次のようになる。

これが、一般線形モデルの一般式となる。ANOVA、単回帰、重回帰、ANCOVAなどはすべてこの一般式で表現が可能である。ここで、Yは測定値の縦ベクトル、Xは説明変数の行列（デザイン行列と呼ばれる）、βは係数の縦ベクトル、εは誤差の縦ベクトル、σは誤差の標準偏差（＝標準誤差）、Iは単位行列である。なおここでは、σは0.55と推定されている。

重要なのが2行目で、これは、「εつまり誤差が平均ゼロ・分散共分散σ2×単位行列の多変量正規分布に従う」ことを示す。言い換えると、「εの9個のどの要素も（つまり、どの測定値の残差も）、互いに独立に、平均ゼロ・分散σ2の正規分布に従う」という意味である。もっと噛み砕いて言えば、「どの測定値の誤差もゼロからσ2程度ランダムにばらつく。また、ある測定値の誤差が上ぶれ/下ぶれしたとき、別の測定値の誤差がそれに連動して上ぶれ/下ぶれすることはない」という意味である。

この2行目における "σ2I" は、このモデルの「誤差の分散共分散行列」と呼ばれる（なお、この行列はしばしばRと表記される）。これがどういうものなのかイメージするために、もう少し具体的に書き下したのが以下である。

上の式のうち、σ2のすぐ右にある9×9の正方行列（行数と列数が同じ）が、このモデルの「誤差のcorrelation structure」である。これは、「ある測定値と別の測定値の誤差の総当たり対戦表」のようなものである。ここでいうcorrelationとは、「似たような値になりやすさ（類似性）」と考えてよい。行、列とも、個体1のTime1, Time2, Time3、個体2のTime1,,,、個体3のTime1,,,のように並んでいる。たとえば、この行列の「1行目、3列目」にある値は、「個体1Time1の測定値（0.86）と個体1Time3の測定値（1.43）の誤差間の類似性（の期待値）」を示す値である。当然、総当たり対戦表なので、「3行目、1列目」にある値も「個体1Time1と個体1Time3の測定値の誤差間の類似性」であり、その他の場所の値もすべて同様なので、この行列は対角線に沿って左右対称の行列である。このように、「誤差のcorrelation structure」は、「総当たり対戦表のような、データ数×データ数かつ対角線に沿って線対称の正方行列」で表される。

さて、このモデルの誤差のcorrelation structureでは、「ある測定値と別の測定値の誤差の類似性は、どのペアであってもゼロ」になっている。これが、ANOVAや線形回帰を含む一般線形モデルが仮定する誤差のcorrelation structureであり、「すべての測定値は互いに独立である」という一般線形モデルの大前提を式で表したものとなっている。

このように、誤差のcorrelation structureは、測定値間の類似性に関する統計モデルの前提を式で表すものである。ということは、誤差のcorrelatoin structureを別の形（値）にすることで、「すべての測定値は互いに独立である」という前提を緩めることができると考えられる。この発想を実現するためのアプローチの一つが、ランダム効果の導入である。

2. ランダム切片モデルのcorrelation structure

2.1. ランダム切片モデル（誤差：正規分布）の作成

先ほどのデータをもう一度確認してみよう。

>dat

     Weight Time Larval_ID

1 0.8625486    1         1

2 1.2923105    2         1

3 1.4263123    3         1

4 1.9185645    1         2

5 2.3913196    2         2

6 2.1835319    3         2

7 1.8324626    1         3

8 1.3362665    2         3

9 1.5063424    3         3

同一個体からの反復測定デザインであるため、同一個体内のデータは、別個体同士のデータに比べて類似する可能性が高いことが考えられる。実際、彼らの体重変化を個体を識別してプロットすると、以下のように、個体1, 3, 2の順に重い傾向にあることがわかる。

ということで、次はセオリー通り個体をランダム効果（ランダム切片）として導入したモデルを作成する。説明のしやすさ上、今回はパッケージ'nlme'を用いたモデル作成を行う。パッケージnlmeでは、lme関数を用いてlinear mixed modelを作成することができる。引数randomを”random=~1|Larval_ID”とすれば、幼虫個体をランダム切片に指定したモデルを作成することができる。結果は以下のようになる。

係数の値は、一般線形モデルの時と同じく、

• Intercept: 1.54

• Time2: 0.14

• Time3: 0.17

となっている。

また、ランダム効果は”標準偏差”で表されており、Interceptが0.46、Residualが0.30となっている。ちなみに、0.462 + 0.302 = 0.552 である。要するにこれは、一般線形モデルの時に推定された誤差分散（誤差の標準偏差の2乗）を2つに分解したものだと解釈できる。

2.2. ランダム切片モデルの行列表示

次に、一般線形モデルのときと同じように、ランダム切片モデルの行列表示を試みる。

ここで、Y, X, β, εはすべて一般線形モデルの時と全く同じであるが、新たにZとbという行列またはベクトルが追加されている。Zは「ランダム効果のデザイン行列」であり、左列から順に、個体1, 2, 3を表すダミーコードとなっている。さらにbは「ランダム効果の係数ベクトル」であり、上から順に個体1, 2, 3の時に切片が全平均から上下する値の大きさを表している。たとえば、一番上の測定値であれば、0.86 = (1×1.54 + 0×0.14 + 0×0.17) + (1×b1+ 0×b2+ 0×b3) + (-0.29) と計算される。

さらに、この行列表示を一般式で表記してみる。

これが、ランダム切片モデルの一般式となる。ここで、Xβをfixed part、Zb + εをrandom partと呼ぶ。Random partとはすなわち、値が確率的に決まる、言い換えれば、測定するたびに値が変わる部分である。σbはランダム切片の標準偏差、σは誤差の標準偏差である。ここで、推定されるパラメータはβの各要素のほか、ランダム切片の標準偏差σb、誤差の標準偏差σである。上記のモデルでは、σb は0.46、σは0.30と推定されている（上記summaryを参照）。

さて、"Y ~ N(BB, CC)"のような表記が何を意味するかはすでに説明した。今回新しく登場した2行目は、bは「平均ゼロ・分散共分散σb2Iの多変量正規分布」に従うことを意味し、σb2Iは「ランダム効果の分散共分散行列」である（しばしばGと表記される）。このモデルの場合、ランダム効果の水準数（使った個体の数）は3なので、行列Gは3×3の正方行列となる。3行目は、一般線形モデルと同じで、σ2Iが誤差の分散共分散行列Rである。

最後に、測定値間の類似性を示す「測定値の分散共分散行列V」を考える。一般線形モデルでは、測定値の類似性＝誤差の類似性なので、行列Rについてだけ考えればよかった。一方、ランダム効果入りモデルでは、測定値の類似性は、誤差の類似性にランダム効果を作用させることで導かれる。行列Vの導出はちょっとややこしいので省略するが、結論だけ示すと、以下のようになる（空白部分はすべてゼロである）。

(σ2 + σb2)の右の行列が、「測定値のcorrelation structure」である。一般線形モデルの行列Rのときと同じく、行、列とも、個体1のTime1, Time2, Time3、個体2のTime1,,,、個体3のTime1,,,のように並んでいる。これを見ると、たとえば行列の1行目、2列目の値がゼロでなくなっていることがわかる。これは、「個体1におけるTime1とTime2の測定値が独立ではない」ということを表している。Correlation structureが総当たり対戦表だということを思い出せば、「測定値のcorrelationがゼロでない（類似している）のは、同一個体からの測定値ペア」であることがわかる。そのcorrelationは、どれも一律で" σb2 / (σ2 + σb2) "である。一方で、別個体同士の測定値のcorrelationはどのペアであってもゼロである。

補足的だが、もう一つ覚えておくとよいのは、行列Vの表記の省略方法である。上記行列Vは、3×3（3時点分）正方行列が対角成分に3つ（3個体分）並んでいるものと見ることができる。これは、同一個体からの測定値の総当たり対戦表である。この3×3正方行列を、Vi（iは個体のID）と定義する。V1=V2=V3、かつ、個体間のcorrelationはゼロなので、測定値のcorrelation structureを考えたいときは、V全体ではなくViだけ見ればよいことになる。

このように、ランダム効果を導入したモデルは、誤差のcorrelation structureにランダム効果を作用させることで、測定値間が独立だという前提を緩めたものだということができる。ランダム切片モデルの場合、ランダム効果の同じ水準（この例では、同じ幼虫個体）から測定された値はすべて同程度に類似している。

2.3. Rによるcorrelation structureの表示

lme関数で作成したランダム効果モデルでは、getVarCov関数によって、これらの分散共分散行列G, R, Vを自動で表示させることができる。（これがlme4パッケージではできないので、nlmeパッケージを使ったのだ）

1つ目の出力が、行列Gを表している。type="random.effects"と指定することで出力される。σb は0.46と推定されているので、表示されるのは0.462 = 0.21となる。なお、この行列はこの値0.21が対角成分に並び、他がすべてゼロの正方行列なので、他の部分はすべて省略され、0.21だけが表示されている。

2つ目の出力が、行列Rである。type="conditional"と指定することで出力される。これは、ランダム効果σは0.30と推定されているので、0.302 = 0.09が対角成分に並び、他がすべてゼロの正方行列が表示されている。ただし、行列Rの全体像はVと同じように3×3の正方行列（Ri）が対角線上に3つ並ぶような行列になる。ここで表示されているのは、そのRiである。

3つ目の出力が、行列Vである。type="marginal"と指定することで出力される。これも、実際に表示されているのは行列Viである。上でも示したViの形を計算すると、対角成分はたしかに" (σ2 + σb2) × 1 = σ2 + σb2 = 0.462 + 0.302 = 0.30 "となり、それ以外の成分は" (σ2 + σb2) × {σb2 / (σ2 + σb2)} = σb2 = 0.462 = 0.21"となることが確かめられる。

3. 誤差のcorrelation structureを直接指定する

3.1. 誤差のcorrelation structureに"compound symmetry"を指定

さて、ここまではセオリー通りのlinear mixed modelのやり方であった。ランダム切片モデルでは、誤差の共分散行列Rにランダム効果を作用させることで、測定値のcorrelation structureを対角成分以外に" σb2 / (σ2 + σb2) "が並ぶように変化させていた。このように、対角成分以外にすべて同一の値が並ぶcorrelation structureを、「compound symmetry」と呼ぶ。

しかし、もし行列Rを直接いじることができれば、ランダム効果など使うことなく、測定値のcorrelation structureを「compound symmetry」の形にできるのでは？とも考えられる。つまり、以下のようなモデルである。ここで、ρは同一幼虫個体から得られたデータの誤差の類似性を示す相関係数である。

このモデルにはランダム効果がないので、誤差の分散共分散行列Rはそのまま測定値の分散共分散行列Vとなる。パッケージnlmeには、行列Rを直接いじるための引数"correlation"が実装されている。以下に例を示す。

このモデルはランダム効果を含まないが、ひとつ前のセクションで用いたlme関数はランダム効果がないと使えないので、ここでは代わりにgls関数を使っている。引数correlationにcorCompSymm(form=~1|Larval_ID)を指定することで、同一幼虫個体について、correlation structureがcompound symmetryであることを指定している。

面白いことに、Fixed effectsの係数推定値、SE、t値が、ランダム切片モデルの値とまったく同じになっている。さらに、AIC、BIC、log likelihoodまですべて同値である。一方、residualは一般線形モデルのものと同値となっている。（t値が同一なのにP値が異なっているのは、ランダム切片モデルが自由度に4を取っているのに対して、correlation: compound symmetryのモデルが自由度に6を取っていることによる。Random partをいじったモデルの自由度をどう取るかには議論があるようだ。）

新たに登場したのがcorrelation structureのセクション中の"Rho"である。これは上記の行列Rにおけるρのことで、0.70と推定されている。したがって、誤差のcorrelation structureは、対角成分に1、それ以外に0.70が並ぶ行列になる。なお、このモデルでは、推定されるパラメータはβの各要素のほか、誤差の標準偏差σ、誤差のcorrelation ρとなる。

さらに、σは0.55と推定されているので、R = Vの対角成分以外の値は、σ2 × ρ = 0.552 × 0.70 = 0.21となるはずである。では、これをgetVarCov関数で確かめてみよう。

3.2. ランダム切片モデルと"correlation=compound symmetry"モデルの行列Vの比較

gls関数で作成したモデルも、getVarCov関数でcorrelation structureを取り出すことができる（ただし、この場合は行列Vしか出力できない）。

狙い通り、まったく同じ行列になっている。

このように、ランダム切片モデルと、行列Rのcorrelation structureにcompound symmetryを指定したモデルは、実質的に全く同じ結論を導くモデルだということができる。

行列Rのcorrelation structureを通して作用する効果のことを、SASの用語では”R-side effects”と呼ぶ。一方、ランダム切片のような通常のランダム効果を通して作用する効果は”G-side effects”と呼ばれる。

4. もっと柔軟なモデリング

4.1. AR-1 correlation structureとランダム切片を混在させる

もう一度、今回のデータを見直してみよう。

>dat

     Weight Time Larval_ID

1 0.8625486    1         1

2 1.2923105    2         1

3 1.4263123    3         1

4 1.9185645    1         2

5 2.3913196    2         2

6 2.1835319    3         2

7 1.8324626    1         3

8 1.3362665    2         3

9 1.5063424    3         3

今回のデータは単に同一個体から複数回データを取っているだけではなく、「経時的に測定を行っている」データである。すると、以下のようなことが考えられる。

• Time1の測定値とTime2の測定値は、測定した時間が近いので、類似する可能性が高い

• Time1とTime3の測定値は、相対的に測定した時間が遠いので、類似性が低くなる

• 「同一個体からの測定値は別個体からの測定値よりも類似している」という効果以外に、「時間的な距離が遠くなるほど類似性が低下する」という効果も考慮する必要がある

「同一個体からの測定値は別個体からの測定値よりも類似している」という効果は、個体をランダム切片に指定することで導入可能である。

一方、「時間的な距離が遠くなるほど類似性が低下する」という効果は、以下のような誤差のcorrelation structureを指定することで考慮することができる。

ここでは、誤差のcorrelation structureを幼虫個体ごとに指定することを仮定しているため、行列Riを表示する。Compound symmetryの場合と同じように、ρは値の類似を示す相関係数である。1行目に注目すれば、左から順に、1, ρ, ρ2と並んでいる。このうち、2つ目はTime1と2のcorrelation、3つ目はTime1と3のcorrelationであり、ρは多くの場合0 < ρ < 1なので、「時間的な距離が遠くなるほど、類似性が低下する」という効果が表現されている。（ひとつ前の時点に影響されて現時点の値が決まる、のようなプロセスを考慮するものというよりは、単に距離によって相関が低下する、というパターンの表現として解釈するのが適切だと思う）

このように、行列上の距離が遠くなるごとに相関係数ρを一つずつ掛けていくようなcorrelation structureを「AR-1」と呼ぶ。このcorrelation structureは、一定間隔の時間的あるいは空間的な位置関係がある時に適用可能である。

関数lmeでは、ランダム効果（G-side effects）と誤差のcorrelation structure（R-side effects）を同時に指定することが可能である。個体をランダム切片に、個体ごとの誤差のcorrelationをAR-1に指定したモデルは以下のようになる。

このモデルの場合、行列Viは以下のようになる。（そろそろ相関行列（correlation structure）で表すのがしんどくなってきたので、分散共分散行列です。行列から (σ2 + σb2) をくくり出せば相関行列になります）

これにより、「同一個体からの測定値間では値が類似する」という効果（σb2）と、「時間的な距離が離れるほど測定値が類似しなくなる」という効果（σ2 ρm）の2つが導入されていることが読み取れる。

それでは、幼虫個体をランダム切片に、誤差のcorrelation structureをAR-1に、同時に指定してみよう。

こちらが、ランダム切片に幼虫個体を、誤差のcorrelation structureにAR-1を指定したモデルの出力となる。誤差のcorrelation structureを指定するには、引数correlation=corAR1(form=~1|Larval_ID)とすればよい。

Correlation structureがAR(1)と表記されるようになり、推定された相関係数は0.44となっている（一般的な教科書では相関係数をρ(rho)、共分散をφ(phi)と表記することが多いが、RのAR-1の出力ではなぜか相関係数なのにphiと表記されている）。このモデルでは、推定されるパラメータはβの各要素のほか、ランダム切片の標準偏差σb、誤差の標準偏差σ、誤差のcorrelation ρの合計6つとなる。

Fixed effectsの係数はここまでのすべてのモデルと同一だが、係数のSEやAICが変わっていることに注意しよう。ランダム切片のみ、またはcompound symmetryのモデルと比べると、Time2のSEは減少しているが、Time3のSEはむしろ増加している。また、AICもやや増加しているため、AR-1 structureの導入はしない方がデータへの当てはまりはよさそうである（経験的には、時点の数が多い方がAR-1的なstructureは当てはまりやすくなると思う）。ついでに、誤差の標準偏差もランダム切片のみのモデルと比べると少し増大している。

4.2. ランダム切片+"correlation=AR1"モデルのG、R、V matrix

最後に、ここまで見てきたstructureの指定が成功しているか、getVarCov関数で確認してみよう。

行列Ri（type="conditional"）、行列Vi（type="marginal"）ともに、行列上の距離が離れるほど値が小さくなっていることが読み取れる。

行列Riの対角成分以外は σ2 ρm なので、 0.392 × 0.441 = 0.07、 0.392 × 0.442 = 0.03であり、行列Viの対角成分以外は σ2 ρm + σb2 なので、 0.392 × 0.441 + 0.382 = 0.21、0.392 × 0.442 + 0.382 = 0.17となる。よって、計算通りにcorrelation structureが導入されたことが確かめられた。

5. 結論のようなもの

ここまで、ランダム効果とcorrelation structureの関係について解説してきた。ランダム切片は、測定値のcorrelation structureをcompound symmetry型に変化させるものであること、誤差のcorrelation structureの指定（R-side effects）を併用することで、ランダム効果（G-side effects）だけを用いるよりも柔軟なモデリングが可能であることなどが分かったと思う。

近年、Rでのmixed modelingのデファクト・スタンダードとなっているパッケージlme4やglmmTMBでは、このcorrelation structure（またはcovariance structure）を可視化することができないため、mixed modelingの文脈におけるcorrelation structureは非常に学習しづらい概念になっている。しかし、理論は本稿で用いたパッケージnlmeと同一である。よって、lme4やglmmTMBの使用で迷ったら、とりあえずnlmeでcorrelation structureを確かめてみるのが良いと思う。Correlation structureを意識しながらモデリングすることで、mixed modelingの間違いはかなり少なくなるはずである。

'21 1/25