反復測定データのモデリング：(5)一般化線形混合モデル

はじめに

前回は、正規分布に従うと想定できないような反復測定データの分析を行うための有力な方法である一般化推定方程式(GEE)を紹介しました。この方法は一般化線形モデル(GLM)と同様のモデリングを行うのですが、目的変数に特定の確率分布を想定せずにパラメータの推定を行うのが特徴です。

今回は、正規分布を仮定せずに分析を行うもう一つの代表的な方法である、一般化線形混合モデル(Generalized linear mixed model; GLMM)について見ていこうと思います。
GLMMはGLMと混合モデルとを組み合わせたモデルであり、GLMを直接拡張したモデルであると言えます。

f:id:mstour:20210323212456j:plain

GLMMもGEEと同じく、整数値のデータ、二値のデータなど、さまざまなデータに対して適用することができます。まずは典型的な使用例を紹介し、続いて一般的な定義を確認していきたいと思います。

GLMMの例

GEEについて紹介した前回の記事でも説明に用いた、てんかんの発作回数に関する反復測定データを分析することを考えてみましょう*1。
この例では、新治療か既存治療かのいずれかを受けた各個人について、複数の期間に分けててんかんの発作回数を記録し、発作の発現率（週あたりの発作回数）について新治療と既存治療とでどの程度の差があるかを評価することが分析の目的です*2。
mstour.hatenablog.com

個人 $i$ の $j$ 回目の観察期における発作回数を $y_{ij}$ とし、これらはそれぞれパラメータ $\mu_{ij}$ のPoisson分布に従うと想定します。Poisson分布に従う確率変数の期待値はパラメータと等しいので、個人 $i$ の期間 $j$ での発作回数の平均（観測された値の単純な平均ではなく、目には見えないが存在するであろうと想定する「真の値」）が $\mu_{ij}$ であると考えていることになります。
また、 $x_i$ を個人 $i$ が受けた治療を示す変数（新治療の場合1、既存治療の場合0）、 $t_{ij}$ を個人 $i$ の観察期 $j$ の長さ（単位は週とします）とします。
このとき、GEEでは、GLMと同様の以下のようなモデル

$\displaystyle \log \mu_{ij} = \beta_0 + \beta_1 x_i + \log t_{ij}$

を考え、モデルのパラメータを推定する際に目的変数 $y_{ij}$ の間に相関があることを考慮することになっていました。

一方で、GLMMでは目的変数どうしに相関があることをモデル式そのものに反映させます。これは、以前に紹介した（正規分布の場合の）混合モデルと全く同じ考え方です。
mstour.hatenablog.com

つまり、同じ個人の測定値どうしは個人差の影響によってある程度似通ったものになると考えるのが自然であるため、個人差を表すランダム効果をモデル式に含めます。
個人差は説明変数と同じく観測値（目的変数）へ影響を与えるものと考え、例えば

$\displaystyle \log \mu_{ij} = \beta_0 + \beta_1 x_i + r_i + \log t_{ij}$

のようなモデルを想定します。 $r_i$ が個人 $i$ に由来する特徴、つまり個人差を表すランダム効果です。
上記のモデルは、個人 $i$ の各観察期 $j = 1, \cdots, T$ それぞれにおける発作回数の真の値 $\mu_{i1}, \cdots, \mu_{iT}$ に共通する部分 $r_i$ があることを想定しており、観測された値 $y_{i1}, \cdots, y_{iT}$ どうしが似ていることを直接表現しています。

f:id:mstour:20210323212512j:plain

見方を変えると、通常のGLMでは説明変数 $x_i$ による影響以外は誤差的なばらつきとして扱われる一方、ランダム効果 $r_i$ を導入することによって、誤差的なばらつきから個人差による影響を分離することができるため、より精度の高い推定が可能になると言えます（この点は、通常の混合モデルでも全く同様です）。

GEEと異なり、目的変数の確率分布としてPoisson分布を明示的に仮定しているため、関心のあるパラメータ $\beta_1$ は最尤法に基づいて推定することができます。

一般的な定義

続いて、GLMMの一般的な定義を見ていきたいと思います。ここでも反復測定データを念頭に置いた説明をすることにしますが、データどうしに相関の生じるその他の場面、例えば同じ地域に住む人どうしの相関を考慮しなければならないようなときでも基本的な考え方は同じです。
表記は前回の記事で紹介したGEEと同様のものを使用することにします。すなわち、目的変数である反復測定データベクトルを $\boldsymbol{y}_i (i = 1, \cdots, N)$ 、それらが従う確率分布の平均パラメータを $E(\boldsymbol{y}_i) = \boldsymbol{\mu}_i$ とします。

通常のGLM（およびGEE）の場合、平均パラメータと説明変数との関係を以下のモデル式

$\displaystyle g(\boldsymbol{\mu}_i) = X_i^T \boldsymbol{\beta}$

で表現しました。ここで行列 $X_i^T$ は、各時点の反復測定データに対応する説明変数ベクトル $\boldsymbol{x}_{ij}^T$ を縦に並べたもの*3です。

GLMMは、このモデルにランダム効果を加えた

$\displaystyle g(\boldsymbol{\mu}_i) = X_i^T \boldsymbol{\beta} + Z_i^T \boldsymbol{u}$

という形のモデル式で表現されます。
$\boldsymbol{u}$ はモデルに含まれる全てのランダム効果からなるベクトルで、 $Z_i$ は「ランダム効果が設定される説明変数」のベクトルを縦に並べたもの*4です。
混合モデルに関する以前の記事でも触れましたが、モデルの切片の部分にランダム効果を設定する*5だけではなく、説明変数の効果に関してもランダム効果を付加することができます*6。 $Z_i$ にはランダム効果が設定された説明変数と、切片に対応する「1」が含まれることになります。

ランダム効果が何らかの確率分布に従うとする点も、正規分布の場合の混合モデルと同様です。ランダム効果（からなるベクトル $\boldsymbol{u}$ ）は（多変量）正規分布に従うものと設定することが一般的です。

このように、GLMMには、目的変数が正規分布の場合の混合モデルと共通する部分が多くありますが、パラメータ（モデル右辺の固定効果 $\boldsymbol{\beta}$ と、ランダム効果・誤差それぞれの分散）の推定において少し異なる部分があります。
正規分布の混合モデルの場合には、例えば固定効果の推定結果*7は

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

のように書き下すことができましたが、GLMMの推定にはさまざまな難点があることが知られています。
そのためGLMMでは、尤度関数を近似することで推定計算を行う方法がいくつか提案されています。

まとめと参考資料

一般化線形混合モデル(GLMM)は、正規分布を想定するのに無理があるような反復測定データの分析に有効な方法の一つです。
前回紹介したGEEと基本的には同様の方法と言えますが、GLMMの説明変数の効果（固定効果）は「同じ人（同じランダム効果を持つ人）が、もし説明変数の値が違っていた場合にどの程度の影響があるか」を示しており*8、このような条件をおかないGEE*9とは多少意味合いが異なるところがあります*10。

少し理論的な内容が続いたため、本シリーズはいったんここで区切り、次回は別の話題を考えたいと思います。

本記事の作成にあたっては以下の資料を参考にしました。
[1] SAS/STAT User's Guide, The GENMOD Procedure.
https://documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=9.4_3.4&docsetId=statug&docsetTarget=statug_genmod_examples07.htm&locale=en
[2] SAS/STAT User's Guide, The GLIMMIX Procedure.
https://documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=9.4_3.3&docsetId=statug&docsetTarget=statug_glimmix_details07.htm&locale=en
[3] RANDOM COEFFICIENT POISSON MODELS | R FAQ, UCLA Statistical Consulting.
https://stats.idre.ucla.edu/r/faq/random-coefficient-poisson-models/
[4] Brown and Prescott(2015), "Applied Mixed Models in Medicine", Wiley.
[5] Gory et al. (2020), "A class of generalized linear mixed models adjusted for marginal interpretability", Statistics in Medicine. 2020;1–14.

*1:データの詳細は、参考資料[1]をご参照ください。

*2:前回も述べた通り、Poisson分布を用いた回帰では、新治療を受けたグループの発現率が既存治療を受けたグループの何倍になるか（率比(Rate ratio)）により結果を要約することが行われます。例えば率比が0.5であれば、新治療は既存治療よりも発作の発現を半減させるという解釈になります。

*3:すなわち、固定効果に関するデザイン行列です。

*4:ランダム効果に関するデザイン行列です。

*5:目的変数の値の平均的な水準に個人差があることを考えていることになります。

*6:介入の効果、時間による変化、リスク因子の影響など、説明変数が目的変数へ与える影響の大きさに個人差があることを表現しています。

*7:正確には、値が確率的に変動する「推定量」のことを指しています。得られたデータから計算される値である「推定値」とは厳密には異なります。

*8:このようなモデルは「条件付きモデル(Conditional model)」や「Subject-specific model」と呼ばれます。

*9:GEEは、個人の効果を考慮せず集団全体での説明変数の効果を考えます。こちらは「周辺モデル(Marginal model)」や「Population-average model」と呼ばれます。

*10:実務上でそこまで厳格に使い分けられているかは、よく分かりません。ただし、GLMMとGEEとでは一般に推定結果は異なります。なお正規分布の場合には、条件付きモデルである混合モデルと、周辺モデルであるCovariance pattern modelとの結果は一致するため、厳格な使い分けが議論になることはあまりないのではないかと思います。