反復測定データのモデリング:(4)一般化推定方程式

はじめに

当ブログではこれまで反復測定データの解析方法のいくつかを紹介してきたが、今のところ目的変数の確率分布に正規分布を想定するものにとどまっていた。
当然、現実の場面ではそうではないデータに対峙することも多く、例えば二値や整数値で記録される反復測定データを扱わなければならない場合もある。

今回は、正規分布以外の反復測定データを解析するための方法として知られている「一般化推定方程式(Generalized estimating equation; GEE)」について考えてみたい。以下、広く用いられている略称であるGEEという表現で統一する。
この方法は、通常の線形モデルとは異なりデータの背後にある確率分布を明示的に特定しないのが大きな特徴である。したがって、正規分布からかけ離れた分布構造を持つデータにも適用可能となるのだが、確率分布を使わないことによるデメリットも存在する。

まずは、GEEの概要を見ていこう。


概要(具体例を交えて)

目的変数が反復測定データではなく1個人につき1個の場合、正規分布以外の確率分布を想定したモデリングを行うには一般化線形モデルを用いることが多い。
例えば、新しい治療が既存治療よりも「てんかん発作」の1週間あたりの発現回数(回/週)を抑えられるかどうかを検討したいという場合、治療開始後の観察期間における発現回数がPoisson分布に従うと想定した一般化線形モデル(Poisson回帰と呼ぶことが多い)を用いることで、発作の発現率(単位時間あたりの発現回数。ここでの例では、週あたりの発現回数を評価したいと考えている)に治療が与える影響を評価することができる。

f:id:mstour:20210312205441j:plain

治療開始から一定の期間にわたって発作発現回数を記録したものとし、その期間の個人 iの発現回数を y_iとする。Poisson回帰によって解析する場合、 y_iはそれぞれパラメータ \mu_iを持つPoisson分布に従うと想定し、 \mu_iについて次のようなモデル

 \displaystyle
\log \mu_i = \beta_0 + \beta_1 x_i + \log t_i

が考えられる。

f:id:mstour:20210312205510j:plain

ここで x_iは治療を表す変数で、新治療を受けた場合1、既存治療を受けた場合0をとる*1。また t_iは個人 iの発作の観察期間(週)とする*2

このモデルにて推定された \beta_1の値の指数をとると、新治療を行った時の「週あたりの発現回数」が既存治療の何倍になるかが算出できる*3

さらに、それぞれの個人のてんかん発作の回数が1つの期間だけではなく、複数の期間に分けて記録されており(つまり、反復測定データになっている)、それらの情報を全て用いて週あたりの発作発現回数を評価したいものとしよう*4
この場合、同じ個人から得られたデータどうしには相関が生じるため、全てのデータを独立とみなしたPoisson回帰モデルをそのまま当てはめるのは適切ではないだろう*5。そこで、一般化線形モデルの理論を踏まえつつ互いに相関のあるデータをモデリングできる方法であるGEEが有用となる。

f:id:mstour:20210312205536j:plain

GEEにおいては、基本的なモデル式は一般化線形モデルと同様のものが用いられる。Poisson回帰を行うここでの事例では、個人 i j回目の測定データ( j回目の観察期間における発作回数) y_{ij}がパラメータ \mu_{ij}を持つPoisson分布に従うと想定し*6、例えば一般化線形モデルの場合と同様に

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

というモデルを考える。
しかしながら、既存治療に対する新治療の効果(週あたり発作回数が何倍になるか)などの未知パラメータの推定は、目的変数に相関があることを考慮した上で行われる。

モデル式のパラメータ(上記の例では \beta_0, \beta_1)をGEEによって推定した結果は、一般化線形モデルの最尤法による推定と同様に真の値に対して一致性を持つ*7。また推定値の分布はサンプルサイズが大きくなるにつれて正規分布に近づくので、このことを利用して関心あるパラメータ(例えば、既存治療に対する新治療の効果を示す \beta_1)についての信頼区間を計算することなどができる。

当然、Poisson分布だけでなくその他の分布が適切と考えられるような場合でもGEEを用いることができる。またSAS、R、SPSSなど広く利用されている統計ソフトにはたいてい実装されているようであり、汎用性の高い方法と言えるだろう。

一方で、GEEは最尤法による推定を行わないため、混合モデルとは異なり欠測メカニズムがMARの場合に推定結果が真の値に近づく(一致推定量になる)ことは保証されていないという難点がある。
欠測メカニズムがMARである場合にも妥当な方法として、一つは一般化線形モデルにランダム効果を導入する「一般化線形混合モデル(Generalized linear mixed model; GLMM)」がある。GLMMは最尤法による推定を行うので、MARの場合にも真の値に近づくような推定結果を得ることができる。
他には、GEEを直接拡張する方法として、観測確率の推定値の逆数で重みづけを行う「IPW推定方程式」と呼ばれる方法を用いることでMARの場合にも妥当な推測ができることが知られている。


数学的詳細

GEEの考え方自体はシンプルで、正規分布の想定が適切とは限らない場合に、互いに相関のあるデータも全て利用することでより精度の高い推定を行うための方法であると理解できる。
しかしながら、数学的な詳細は(少なくとも筆者には)直感的になかなか理解しにくい部分がある。理論的な側面を少しだけ追ってみよう。

まずは通常の一般化線形モデルから話を進める。
関心のある目的変数を y_i (i = 1, \cdots, N)、その期待値を E(y_i) = \mu_iとする。先のPoisson回帰の例では \mu_i y_iが従うと想定するPoisson分布のパラメータである。
これに対して、説明変数ベクトル \boldsymbol{x}_i(切片に対応する1を要素に含む)を用いた一般化線形モデル

 \displaystyle
g(\mu_i) = \boldsymbol{x}_i^T \boldsymbol{\beta}

を当てはめる状況を考えよう。ここで \boldsymbol{\beta} = (\beta_0, \beta_1, \cdots, \beta_P)^Tは未知のモデルパラメータとする。一般化線形モデルについては以前の記事もよければご覧いただきたい。

パラメータの推定を行うには、想定する確率分布から作られる「尤度(Likelihood)」を最大化するようなパラメータの値を求める必要がある。実際には、尤度の対数をとった「対数尤度(Log likelihood)」が最大になるようなパラメータの値を計算する。
対数尤度を \ellとすると、パラメータ \beta_pの推定値は \ell \beta_p偏微分した結果を0と等しいとした方程式を解くことで得られる:

 \displaystyle
\frac{\partial \ell}{\partial \beta_p} = 0

左辺は \beta_pについての「スコア統計量(Score statistics)」といい、以降では U_pと表す。またスコア統計量が0に等しいとおいた方程式は「スコア方程式(Score equation)」と呼ばれ、以下のように表されることが知られている。

 \displaystyle
U_p = \sum_{i=1}^{N} \frac{(y_i - \mu_i)}{Var(y_i)} \frac{\partial \mu_i}{\partial \beta_p} = 0 \hspace{10pt} (p = 0, 1, \cdots, P)

全てのパラメータに関するスコア統計量を縦に並べたものを Uとすると、次のように表せる。

 \displaystyle
U = \sum_{i=1}^{N} \frac{(y_i - \mu_i)}{Var(y_i)} \frac{\partial \mu_i}{\partial \boldsymbol{\beta}} = \boldsymbol{0} \hspace{30pt} (1)


この方程式を \boldsymbol{\beta}について解いた結果がパラメータの真値の一致推定量になることは最尤法の理論から導かれるが、実は目的変数の従う確率分布を具体的に特定しなくても、ある条件を満たせば同様の結果が得られることが示されている。
確率分布を特定しない場合に、上記のスコア方程式は「擬似尤度推定方程式(Quasi-likelihood estimating equation)」と呼ばれている*8

これをさらに一般化して、確率分布を特定せず、さらに互いに相関のあるデータの場合にも拡張したものがGEEである。
GEEを用いる場面を想定して、今度は目的変数を反復測定データのベクトル \boldsymbol{y}_i (i = 1, \cdots, N)とする。目的変数は、平均ベクトル E(\boldsymbol{y}_i) = \boldsymbol{\mu}_i、共分散行列 V_iとなるような何らかの(多変量)確率分布に従うものとする*9。また、平均に関して一般化線形モデル

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

によるモデリングを考える(行列 X_i^Tは、各時点の反復測定データに対応する説明変数ベクトル \boldsymbol{x}_{ij}^Tを縦に並べたもの)。

GEEによるパラメータ推定は、式(1)を多変量に拡張した

 \displaystyle
U = \sum_{i=1}^{N} \left( \frac{\partial \boldsymbol{\mu}_i}{\partial \boldsymbol{\beta}} \right) V_i^{-1} (\boldsymbol{y}_i - \boldsymbol{\mu}_i) = \boldsymbol{0} \hspace{30pt} (2)

をパラメータベクトル \boldsymbol{\beta}について解くことで行う*10
もっとも、実際に解くためには反復計算が必要となる。

データ間に相関があることを考慮するために、式(2)の V_iに、想定されるデータ間の相関の構造を直接指定することができる。やり方は以前に紹介したCovariance pattern modelと同様である。

GEEによるパラメータの推定結果は一致性を持つことが知られているが、例えば分散の推定方法など他にも考えないといけない点がある。長くなってきたので、今度また応用編的な内容として続きを書いていきたい。


まとめ

今回は、正規分布ではモデリングできないような反復測定データにも使用できる方法である一般化推定方程式(GEE)をざっと紹介した。
多くの統計ソフトで実行可能であり、特に整数値をとるカウントデータや二値データの解析に広く活用されているようである。詳しい話はまだ書ききれていないので、また次回へ続くということで・・・


参考文献
[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] Yang et al.(2020), "Transcranial direct current stimulation reduces seizure frequency in patients with refractory focal epilepsy: A randomized, double-blind, sham-controlled, and three-arm parallel multicenter study", Brain Stimul. Jan-Feb 2020;13(1):109-116.
https://www.brainstimjrnl.com/article/S1935-861X(19)30369-9/fulltext
[3] Annette J. Dobson(2008), "一般化線形モデル入門", 共立出版.
[4] John M. Lachin(2020), "医薬データのための統計解析", 共立出版.
[5] 高井啓二・星野崇宏・野間久史(2016), "欠測データの統計科学", 岩波書店.

*1:交絡因子などを含めると話が複雑になるので、治療はランダムに割り付けられているものとする。

*2: \log t_iはオフセット項と呼ばれる。例えば観察期間が長い個人ほど発作の発現回数は多くなる可能性が高いので、全員の観察期間を揃えて評価するために使用される。少し古い記事だがこちらでも紹介した。

*3:この指標は、率(Rate)の比をとったものなので「率比(Rate ratio)」と呼ばれることがある。

*4:ここでの例は文献[1]に記載の古典的な例を参考にしたが、例えば文献[2]など最近でも同様の事例が見つけられる。

*5:文献[1]では推定結果のばらつき(標準誤差)を小さく見積もってしまう例が紹介されているが、文献[3]では逆にばらつきを大きく見積もってしまう例がある。

*6:ただし一般化線形モデルと異なり、確率分布が推定計算の途中で出てくることはない。

*7:つまり、サンプルサイズが大きくなれば真の値に近づいていく。

*8:確率分布を使っていないので尤度に基づいた推定ではないが、擬似的に最尤法と同じ手続きを踏むことからこのような名称になっているかと思われる。

*9:具体的な分布は決めなくてもいいが、平均と分散は何らかの方法で推定できることを前提としている。

*10:総和記号の直後の偏微分は、  \boldsymbol{\mu}_iの要素を \boldsymbol{\beta}の要素で微分したものが並んだ行列になる。