線形混合モデルにおけるKenward-Roger法：(1)概要

【2021/5/10 追記】
Kenward-Roger法による分散推定量の数式について、より正確な記載に改めました。

はじめに

データどうしの相関を考慮したモデルである混合モデルは、反復測定データの解析などで有効であることが広く知られています。
特に、目的変数の分布として多変量正規分布を想定する線形混合モデル*1が用いられる事例は例えば臨床試験のデータ解析の分野では頻繁に見られます。

実際に線形混合モデルを応用する際には、共分散行列の構造（あるいはランダム効果を置く場所や確率分布）などモデルを特徴づけるいくつかの設定を行うのですが、そのうちの一つに「自由度の計算方法」というものがあります。
ここでいう自由度とは、関心のある固定効果（例えば、介入Aの効果）あるいはそれらの「対比」（例えば、ある時点における介入Aと介入Bとの間の差を表現するのに用います。後述します。）の検定を行う際の、検定統計量が従う確率分布*2の自由度を指します。

反復測定データにありがちな欠測が生じた場合でも、混合モデルを用いれば欠測のある個体のデータも利用してパラメータの推定を行うことが可能ですが、自由度の計算は単純ではなくなることが知られています。
自由度の計算方法の選択肢はいくつかありますが、その中の1つ「Kenward-Roger法」の使用が推奨されることが多くあります（例えば参考資料1）。

今回はこのKenward-Roger法について、少し見ていきましょう。
なお、提案者であるKenward and Rogerの原論文（参考資料2）に沿って話を進めていきますが、非常に難しい内容で細かな部分まで理解しきれない*3ので、実用上知っておくといいかなという程度にかいつまんで説明してみたいと思います。

長くなるので、具体的な数値例などを続編として次回また書きたいと思います。

固定効果の分散が過小に見積もられる問題

今回は、以下の線形混合モデル
$\displaystyle \boldsymbol{y} = X \boldsymbol{\beta} + Z \boldsymbol{\gamma} + \boldsymbol{\varepsilon}$

を考えます。 $\boldsymbol{\beta}$ は固定効果のベクトル、 $\boldsymbol{\gamma}$ はランダム効果のベクトルです。
このとき目的変数ベクトル $\boldsymbol{y}$ の確率分布は多変量正規分布 $N(X \boldsymbol{\beta}, \Sigma)$ で表すことができます。

通常、固定効果 $\boldsymbol{\beta}$ や共分散行列 $\Sigma$ に含まれる未知パラメータ*4の推定には制約付き最尤法(Restricted maximum likelihood method; REML)が用いられます。REMLにより推定した共分散行列を $\hat{\Sigma}$ とします*5。以下では話の流れによって「共分散行列」を単に「分散」と言う場合があります。

一方、固定効果を求めるために最適化計算をすると推定量として
$\displaystyle (X^T \Sigma^{-1} X)^{-1} X^T \Sigma^{-1} \boldsymbol{y}$

が得られますが、未知のパラメータを含む $\Sigma$ が存在するため、以下のように推定量 $\hat{\Sigma}$ に置き換えたものが固定効果 $\boldsymbol{\beta}$ の推定量として用いられます。
$\displaystyle \hat{\boldsymbol{\beta}} = (X^T \hat{\Sigma}^{-1} X)^{-1} X^T \hat{\Sigma}^{-1} \boldsymbol{y}$

この $\hat{\boldsymbol{\beta}}$ によって真の値 $\boldsymbol{\beta}$ をバイアスなく推定できることが知られています。

また $\hat{\boldsymbol{\beta}}$ の、サンプルサイズを無限に大きくしたときの（=漸近分布の）分散は
$\displaystyle V(\hat{\boldsymbol{\beta}}) \cong (X^T \Sigma^{-1} X)^{-1} \hspace{20pt} (1)$

となりますが、やはり未知の $\Sigma$ が含まれているため、慣習的に $\hat{\Sigma}$ で代用した
$\displaystyle \hat{V}(\hat{\boldsymbol{\beta}}) = (X^T \hat{\Sigma}^{-1} X)^{-1} \hspace{20pt} (2)$

を分散の推定量として用いてきたようです。
しかしながらこの方法は

$\hat{\Sigma}$ のばらつきが考慮されていない（・・・ $\hat{\boldsymbol{\beta}}$ も $V(\hat{\boldsymbol{\beta}})$ も本来は真の分散 $\Sigma$ を使って算出されるところ、同じ実験を繰り返すと毎回結果が変わる $\hat{\Sigma}$ で代用するため、その分ばらつきが大きくなるはずです）
そもそも $(X^T \hat{\Sigma}^{-1} X)^{-1}$ は $(X^T \Sigma^{-1} X)^{-1}$ の不偏な推定量ではない（・・・式(2)では式(1)を正しく推定できていないということでしょう）

ことから、サンプルサイズが小さい場合には $\hat{\boldsymbol{\beta}}$ の分散を小さめに見積もってしまうことが知られていました。

そこで、KenwardとRogerにより、新たな分散推定量が提案されました。

Kenward-Roger法による分散推定量

固定効果推定量 $\hat{\boldsymbol{\beta}}$ の正確な分散（サンプルサイズを無限に大きくしたときに収束する式(1)とは異なる）、以下のように表現できることが知られていました。
$\displaystyle V(\hat{\boldsymbol{\beta}}) = (X^T \Sigma^{-1} X)^{-1} + \Lambda \hspace{20pt} (3)$

前述のようにサンプルサイズが無限に大きくなると分散は右辺第1項に等しくなりますが、逆に言うとサンプルサイズが小さい場合第1項を用いるだけでは $\Lambda$ の分ほど分散を過小評価してしまいます。

Kenward-Roger法は、分散の過小評価を防ぐために、前述の式(2)の推定量
$\displaystyle \hat{V}(\hat{\boldsymbol{\beta}}) = (X^T \hat{\Sigma}^{-1} X)^{-1}$

に、式(3)の $\Lambda$ に対応する修正項を加えたものを推定量とする方法です。

途中の導出は省略します*6が、式(2)の既存の推定量を $\hat{\Phi}$ とすると

$\displaystyle \hat{V}_A(\hat{\boldsymbol{\beta}}) = \hat{\Phi} + 2\hat{\Phi} \left\{ \sum_{i=1}^{r} \sum_{j=1}^{r} W_{ij} (Q_{ij} - P_i \hat{\Phi} P_j - \frac{1}{4}R_{ij}) \right\} \hat{\Phi}$

がKenward-Roger法による分散推定量となります。ただし $r$ は共分散行列に含まれる未知パラメータ数で、
$\displaystyle P_i = X^T \frac{\partial \Sigma^{-1}}{\partial \sigma_i} X, \\ Q_{ij} = X^T \frac{\partial \Sigma^{-1}}{\partial \sigma_i} \Sigma \frac{\partial \Sigma^{-1}}{\partial \sigma_j} X, \\ W_{ij} = V\left[(\hat{\sigma}_1, \cdots, \hat{\sigma}_r) \right]のi行j列目の要素,$

$\displaystyle R_{ij} = X^T \Sigma^{-1} \frac{\partial^2 \Sigma}{\partial \sigma_i \partial \sigma_j} \Sigma^{-1} X$

です*7。ただし、実際の計算では $\Sigma$ は $\hat{\Sigma}$ で代用し、 $\Sigma$ の微分の項は微分した後の結果に各パラメータ $\sigma_i$ の推定量を代入します。

なお共分散構造がパラメータに関して線形である場合、すなわち
$\displaystyle \Sigma = \sum_{i=1}^{r} \sigma_i G_i \hspace{20pt} (G_iは適当なn \times n行列)$
と表せる場合には、2回微分すると0になるため $R_{ij} = 0$ となり、Kenward-Rogerの分散推定量は次のようになります。

$\displaystyle \hat{V}_A(\hat{\boldsymbol{\beta}}) = \hat{\Phi} + 2\hat{\Phi} \left\{ \sum_{i=1}^{r} \sum_{j=1}^{r} W_{ij} (Q_{ij} - P_i \hat{\Phi} P_j) \right\} \hat{\Phi}$

例えば、共分散構造がCompound symmetryやUnstructuredの場合にはこれに該当しますが、AR(1)のような構造だと相関パラメータの2乗以上の要素が共分散行列に含まれるため、より一般的な $R_{ij}$ を含む推定量の計算が必要になります（SASではこちらがデフォルトのようです）。

f:id:mstour:20210505094228j:plain

統計ソフトウェアの説明などではKenward-Roger法は「自由度計算方法」として記載されていることが多くみられますが、方法が提案された背景としてまずは分散を過小評価する問題があったようです。

固定効果に関する検定の自由度の計算

モデルの固定効果に関して検定を行う場合、対比(Contrast)を使って考えることが多くあります。
対比とは、固定効果パラメータの線形結合（各パラメータに何らかの係数を掛けて足し算したもの。ただし係数の合計が0になるように制約をつける）のことをいいます。

例えば治療A、B、Cの効果を比較する臨床試験を行い、治療AとBとの比較に関心があるものとします。
統計モデルとして全体平均 $\mu$ と、治療A〜Cのそれぞれのアウトカムへの影響を示す固定効果パラメータ $\beta_A, \beta_B, \beta_C$ の4個の固定効果パラメータを含む線形モデルを用いたものとすると、治療AとBの比較 $\beta_A - \beta_B$ は
$\displaystyle L^T \boldsymbol{\beta} = (0, 1, -1, 0) \left( \begin{array}{c} \mu \\ \beta_A \\ \beta_B \\ \beta_C \end{array} \right) = \beta_A - \beta_B$

のように表現できます。 $L$ は対比の係数からなる行列で、上記のように関心ある比較が1つの場合には単にベクトルになります。
上のように対比を定義すると、治療AとBの間に差がないという帰無仮説 $H_0: \beta_A - \beta_B = 0$ は $H_0: L^T \boldsymbol{\beta} = 0$ と表せます*8。

検定には、Wald 検定と呼ばれる方法が用いられます。
対比を用いて表現すると、帰無仮説 $H_0: L^T \boldsymbol{\beta} = 0$ を検定するための統計量は
$\displaystyle F = \frac{(L^T \hat{\boldsymbol{\beta}})^T (L^T \hat{V}(\hat{\boldsymbol{\beta}}) L)^{-1} (L^T \hat{\boldsymbol{\beta}})}{r}$

であり、これが自由度 $r, \nu$ のF分布に近似的に従うことを利用します。ただし $r$ は検定の対象となる（仮説として設定する）対比の数であり、これが1つの場合には以下のt統計量
$\displaystyle t = \frac{L^T \hat{\boldsymbol{\beta}}}{\sqrt{L^T \hat{V}(\hat{\boldsymbol{\beta}}) L}}$

が近似的に自由度 $\nu$ のt分布に従うことを用いて検定することができます*9。
近似分布であるF分布（あるいはt分布）の自由度 $\nu$ は何らかの方法で推定する必要があり、古くからSatterthwaiteにより提案された方法が知られていました。

前置きが長くなりましたが、その自由度の推定方法がKenward-Roger法のもう一つの内容です。
ここまで書いておいて何ですが、自由度推定方法は分散推定量の導出以上に僕の理解できるレベルをはるかに超えているので・・・要点だけ書くと

Satterthwaiteの方法で計算に使用する固定効果の分散推定量として、Kenward-Rogerの提案した過小評価を防ぐ推定量を用いる（参考資料4）
複数の対比の場合には非常に複雑な問題になるが、対比の数が1つの時には、Satterthwaiteの方法の単純な修正である（参考資料2）

ということのようです。なお、検定統計量Fやtの計算にもKenward-Roger提案の分散推定量を使用します。

この点は次回、実際にソフトウェアで計算してどのような違いがあるかを見ていくことにします。
原論文を見ると、Kenward-Roger以前の方法だとType1 errorの確率が設定した有意水準を大きく超えるケースがあるようです。

おわりに

Kenward-Roger法について、原論文を片手に概要を紹介していきました。
非常にざっくりとした解説になりましたが、ポイントは以下の2点にまとめられるかと思います。

サンプルサイズが小さいときには固定効果の分散を過小推定してしまう危険があるので、それを防ぐための適切な推定量が提案された
固定効果の検定に使用する近似分布の自由度の計算でも、新たに提案した分散推定量を使用することにした

かなり数学的に高度な内容で、今ひとつスッキリと理解できませんが、具体的な計算をしてみるともう少し見えてくることがあるかもしれません。
ということで、次回につづきます・・・。

今回は以下の資料を参考にしました。
[1] 日本製薬工業協会(2016). 欠測のある連続量経時データに対する統計手法について Ver2.0.
http://www.jpma.or.jp/medicine/shinyaku/tiken/allotment/statistics.html
[2] Kenward, M., & Roger, J. (1997). Small Sample Inference for Fixed Effects from Restricted Maximum Likelihood. Biometrics, 53(3), 983-997.
[3] Brown and Prescott(2015). Applied Mixed Models in Medicine. Wiley.
[4] SAS/STAT User's Guide: The MIXED Procedure.
https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_mixed_syntax10.htm#statug.mixed.modelstmt_ddfm
https://documentation.sas.com/doc/en/pgmsascdc/9.4_3.3/statug/statug_mixed_details01.htm#statug.mixed.mixedinference
[5] Kuznetsova, A., Brockhoff, P., & Christensen, R. (2017). lmerTest Package: Tests in Linear Mixed Effects Models. Journal of Statistical Software, 82(13), 1 - 26.
https://www.jstatsoft.org/article/view/v082i13

*1:正規(Normal)やガウス(Gaussian)といった言葉を付け加えた方がより正確ですが、「線形混合モデル」で通じることも多いので、煩雑さを避けてこのように記載することにします。

*2:1つの対比についての検定をする場合にはt分布の自由度ですが、より一般の場合（複数の対比の検定）にはF分布の2番目の自由度（分母の自由度）です。

*3:あくまで僕にとってですが・・・

*4:例えば、共分散行列の構造にCompound Symmetry型を仮定する場合には未知パラメータは「各変数で共通の分散」と「すべての変数のペアで共通の共分散」の2つ、Unstructured型の場合には「変数ごとに異なる分散」と「すべての変数ペアごとに異なる共分散」です。

*5:原論文では、共分散行列に含まれるパラメータを $\sigma$ 、その推定量を $\hat{\sigma}$ とし、共分散行列自体はそれらに依存するという意味で $\Sigma(\sigma)$ 、 $\Sigma(\hat{\sigma})$ と表記しているのですが、ややこしくなるので本文のような表記法とします。