反復測定データのモデリング：(1)Covariance pattern model

はじめに

重回帰分析や一般化線形モデルといった線形モデルでは通常、異なる目的変数どうしは独立であることを仮定している。

例えば、B高校の2年生の期末テストの得点を目的変数としたい場合を考えてみると、各学生の得点どうしは独立である*1と仮定してもだいたい問題はないかと思われる。
一方で、B高校2年生であるCさんの「毎月月末に行われる学力テストの得点」についてはどうだろうか。きっと、Cさんの6月の得点と7月の得点は似たような結果になるだろう。すると、これらの得点をあたかも異なる学生の点数のように独立な値と考えるのは少し無理があると考えられる。
f:id:mstour:20210205211205j:plain

反復測定データとは

この例のように、同じ個人から時間とともに何度も繰り返し測定して得られたデータは「経時測定データ(Longitudinal data)」と呼ばれる。
また別の例として、A県内の各高校で実施された共通学力テストの結果に関心があるものとしよう。すると、入試のいわゆる偏差値が高い高校ほど学力の高い学生が集まっていることは容易に想像できるため、同じ高校の学生どうしの点数は似たようなものになると考えるのが自然である。このようにある集団内の類似した個人から得られたデータは「クラスターデータ(Clustered data)」といわれる。

経時測定データ、クラスターデータはいずれも「反復測定データ(Repeated measures data)」と表現されることが多い。つまり、ある単位（経時測定データの例では「個人」、クラスターデータの例では「特定の集団」）から何度も反復してデータを測定していることを意味しており、そのようなデータどうしには相関が生じる。
したがって、目的変数どうしに相関があることを考慮したモデルを使うことが必要となる。

いくつかの選択肢があるが、ここでは「Covariance pattern model*2」と呼ばれるモデルについて見ていこうと思う。
なおこのモデルは正規分布のデータを解析する文脈で用いられることが多いように思われるため、この記事でも同様に目的変数は正規分布に従うものとする。ただし、反復測定されたデータの相関を考慮に入れる必要があるので、データ1個1個が1変量の正規分布に従うとするのではなく、ある個人（クラスターデータの場合はある集団）の反復測定データ全てが多変量正規分布に従うと考える。

モデルの定式化

以下では、同じ個人から繰り返し測定された経時測定データで考えることにする。
いま、各個人 $i(i = 1, \cdots, N)$ について、それぞれ $n_i$ 回の反復測定がなされているとする（測定回数は個人ごとに違っていてもよい）。個人 $i$ の反復測定データ全体をベクトル

$\displaystyle \boldsymbol{y}_i = (y_{i1}, \cdots, y_{i, n_i})^T$

としてまとめ、さらにそれを全員分縦に並べたものを目的変数ベクトル

$\displaystyle \boldsymbol{y} = (\boldsymbol{y}^T_1, \cdots, \boldsymbol{y}^T_N)^T$

とする。このベクトルの長さは $M = \sum_{i=1}^{N} n_i$ であり、Covariance pattern modelでは次の（M次元）多変量正規分布に従うものと想定する：

$\displaystyle \boldsymbol{y} \sim N(\boldsymbol{\mu}, V)$

$\boldsymbol{\mu}$ は全ての個人の反復測定データ全ての平均を並べたベクトルであり、 $V$ は同じく反復測定データ全体の分散共分散行列を表す。通常は、異なる個人間の測定値は独立と想定することが多いため、以下のようなブロック対角行列を考える：

$\displaystyle V = \left( \begin{array}{cccc} V_1 & O & \ldots & O \\ O & V_2 & \ldots & O \\ \vdots & \vdots & \ddots & \vdots \\ O & O & \ldots & V_N \end{array} \right)$

$V_i(i = 1, \cdots, N)$ は個人 $i$ の測定値だけの（すなわち $\boldsymbol{y}_i$ の）分散共分散行列、 $O$ は全ての要素がゼロの行列で、つまり上記の $V$ は「同じ個人のデータどうしには相関があり」「異なる個人のデータどうしには相関がない」ことを表現している。

さて、重回帰モデルと同様に、目的変数の平均ベクトル $\boldsymbol{\mu}$ を説明変数を用いた線形モデルで表現することを考えよう。P個の説明変数をモデルに含めるものとすると、切片を含めたパラメータのベクトル

$\displaystyle \boldsymbol{\beta} = (\beta_0, \beta_1, \cdots, \beta_P)^T$

と、各個人の反復測定データに関するデザイン行列 $X_i$ を並べた行列

$\displaystyle X = (X_1^T, \cdots, X_N^T)^T$

を用いて平均ベクトルは次のように表される：

$\displaystyle \boldsymbol{\mu} = X \boldsymbol{\beta}$

なおデザイン行列は、1列目には全て1が並び（切片に対応）、2列目〜P+1列目にはそれぞれ1番目〜P番目の説明変数の値が配置される。行方向には、上から順に1回目の測定データ、2回目の測定データ、・・・、に対応する値が並ぶことになる。
簡単な例で説明する。目的変数は「6月、7月、8月それぞれの学力テストの得点」であり、説明変数として「性別を示すダミー変数（男性であれば1）」「6月であることを示すダミー変数（6月であれば1）」「7月であることを示すダミー変数（7月であれば1）」があるものとすると、男子学生 $S$ さんのデザイン行列は

$\displaystyle X_S = \left( \begin{array}{cccc} 1 & 1 & 1 & 0 \\ 1 & 1 & 0 & 1 \\ 1 & 1 & 0 & 0 \end{array} \right)$

となる。

目的変数の確率分布が特定されているので、パラメータの推定には最尤法を用いればよい（正規分布なので最小二乗法でもよい）。つまり、対数尤度をパラメータで微分する（微分したものを0に等しいとする）ことで計算することができる。
結果だけ示すと、推定量 $\hat{\boldsymbol{\beta}}$ は

$\displaystyle \hat{\boldsymbol{\beta}} = (X^T V^{-1} X)^{-1} X^T V^{-1} \boldsymbol{y}$

となる。
この式をみると目的変数ベクトルの分散共分散行列 $V$ が含まれているため、先に何らかの方法で $V$ の推定を行う必要がある。
具体的な推定方法はここでは詳しく述べないことにするけれど、その前の話として、重回帰モデルなどと異なり分散パラメータが1つではないため、まずは $V$ の具体的な形を決めてあげなければならない。分散共分散行列の具体的な形は「共分散構造(Covariance pattern)」と呼ばれる。

共分散構造

共分散構造を決めるとは、つまり、「同じ個人の反復測定データどうしがお互いにどういう関係性にあると考えられるか」をある程度事前に決めておくということである。具体的な推定値はデータから計算するわけだが、その構造はモデルの定義としてあらかじめ決めておかないと話が先に進まない*3。
例えば、6月の成績と7月の成績はほとんど同じになるけれど、6月の成績と12月の成績はほとんど関係がないかもしれない。そういうことが事前の知識から予想できるのであれば、そのことを共分散構造に反映させることでより適切なモデル化ができるだろう。

共分散構造として、以下のものが比較的よく使われる。
なお、実用的には全ての個人で分散共分散行列 $V_i$ は同じと考えることが多いので*4、以下では $V_i$ を代表として表記する。全体の分散共分散行列は同じ $V_i$ を対角線上に並べたものになる。

(1) Unstructured （Generalとも呼ばれる）
これは、全てのデータ間の共分散と全てのデータの分散が異なることを許容する、最も制約の弱い構造である。

$\displaystyle V_i = \left( \begin{array}{ccccc} \sigma^2_1 & \theta_{12} & \theta_{13} & \ldots & \theta_{1, n_i} \\ \theta_{12} & \sigma^2_2 & \theta_{23} & \ldots & \theta_{2, n_i} \\ \theta_{13} & \theta_{23} & \sigma^2_3 & \ldots & \theta_{3, n_i} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \theta_{1, n_i} & \theta_{2, n_i} & \theta_{3, n_i} & \ldots & \sigma^2_{n_i} \end{array} \right)$

したがって、最も柔軟にデータにフィットしたモデル化ができるが、測定回数が増えるとパラメータ数が膨れ上がってしまうため、推定が不可能になる（反復計算を行うので、計算が収束しなくなる）恐れがある。
そのため、計算が収束しなかった場合にはより単純な構造を採用する、といったオプションを前もって決めておくこともよく行われる。

(2) Toeplitz
Unstructuredよりも制約を強くして、全てのデータの分散は等しく、また測定の間隔が同じデータの組（例えば、6月と7月の成績と、8月と9月の成績）の共分散（したがって相関も）は全て同じであるとする。

$\displaystyle V_i = \left( \begin{array}{ccccc} \sigma^2 & \theta_{1} & \theta_{2} & \ldots & \theta_{n_i -1} \\ \theta_{1} & \sigma^2 & \theta_{1} & \ldots & \theta_{n_i -2} \\ \theta_{2} & \theta_{1} & \sigma^2 & \ldots & \theta_{n_i -3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \theta_{n_i -1} & \theta_{n_i -2} & \theta_{n_i -3} & \ldots & \sigma^2 \end{array} \right)$

隣りあう時点のデータどうしよりもお互いに離れたデータどうしの方が関連性が薄くなるというケースはよくあると考えられる。Toeplitz型はUnstructuredほどの柔軟性はないが、そのことを反映できる構造になっている。
f:id:mstour:20210205211228j:plain

(3) First-order autoregressive
この構造もToeplitzと同じくデータの間隔が異なれば関連性も異なると想定しているが、下記のように制約はかなり強くなっている：

$\displaystyle V_i = \left( \begin{array}{ccccc} \sigma^2 & \sigma^2 \rho & \sigma^2 \rho^2 & \ldots & \sigma^2 \rho^{n_i -1} \\ \sigma^2 \rho & \sigma^2 & \sigma^2 \rho & \ldots & \sigma^2 \rho^{n_i -2} \\ \sigma^2 \rho^2 & \sigma^2 \rho & \sigma^2 & \ldots & \sigma^2 \rho^{n_i -3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \sigma^2 \rho^{n_i -1} & \sigma^2 \rho^{n_i -2} & \sigma^2 \rho^{n_i -3} & \ldots & \sigma^2 \end{array} \right)$

$\rho$ は隣りあうデータどうしの相関係数であり、相関係数は-1から1までの範囲の値をとるため、測定の間隔が離れれば離れるほど関連性が薄くなっていくことを表現できる。

(4) Compound symmetry
異なる測定時点のデータどうしの関連性は全て同じであるとする構造であり、非常に強い制約を置いている。

$\displaystyle V_i = \left( \begin{array}{ccccc} \sigma^2 & \theta & \theta & \ldots & \theta \\ \theta & \sigma^2 & \theta & \ldots & \theta \\ \theta & \theta & \sigma^2 & \ldots & \theta \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \theta & \theta & \theta & \ldots & \sigma^2 \end{array} \right)$

過度に複雑な構造よりも推奨されることもあるが、上記(1)〜(3)のようにデータの間隔によって関連性が変わってくる様相は捉えられないので、この構造でよいかどうかは慎重に考えた方がいいかもしれない。

まとめ

同じ個人や、特定の同質な集団から測定された反復測定データは互いに相関が生じる。そのような場合にデータどうしを独立であると仮定する通常の回帰分析は適切とはいえないため、相関があることを考慮したモデルが必要である。
その一つであるCovariance pattern modelは、データの集まりの従う確率分布として多変量正規分布を想定し、データどうしの関連性（分散共分散行列）の構造も含めてモデル化を行う方法である。
これによって、通常の回帰分析と同様に説明変数が目的変数へ与える影響を評価することができる。

今回は以下の文献を参考にした。
[1] Annette J. Dobson(2008), "一般化線形モデル入門", 共立出版.
[2] Brown and Prescott(2015), "Applied Mixed Models in Medicine", Wiley.

*1:以下、統計学的な意味での「独立」を意味するものとする。ざっくりいうと、一方の値が他方の値とは全く関係しないということを指す

*2:適切な日本語訳がなさそうなので、英語での表現を用いることにする

*3:もちろん一つの共分散構造に決めうちする必要はなくて、たくさんの候補の中から最も適切なものをAICなどの基準によって選ぶということも普通に行われている。ただしその場合でも、候補をあらかじめ特定しておかなければ何もできない

*4:全ての個人で異なるとすると、推定すべきパラメータの数が膨大になる