反復測定データのモデリング:(3)ランダム係数モデルについて少し詳しく

はじめに

前回(反復測定データのモデリング:(2)混合モデル)は「混合モデル」の概要と、モデルの切片(目的変数に関する各個人の平均的な水準を表す)だけにランダム効果を設定する「ランダム効果モデル」についてを主に述べた。
今回は、前回説明を先送りにしていた「ランダム係数モデル(Random coefficients model)」について詳しく見ていこう。

より発展的なモデル、例えば目的変数の分布を正規分布以外に広げた「一般化線形混合モデル(Generalized linear mixed model)」や、目的変数を非線形関数によって直接モデル化する「非線形混合効果モデル(Non-linear mixed model)」といった方法も、基本的な考え方は正規分布を仮定するランダム効果モデルやランダム係数モデルと同じなので、ランダム係数モデルを知っておくと分析手法の選択肢が広がるかと思う*1


ランダム係数モデルの概要

ランダム係数モデルとは、切片だけでなく説明変数の係数にもランダム効果を設定したモデルである。

説明のための例として、N人の個人からそれぞれT回測定した反復測定データ y_{ij}(i = 1, \cdots, N; j = 1, \cdots, T)を、時間 t_{ij}を説明変数としたモデルで表現する場合を考える。つまり、反復測定データと時間との関係性を直線で表現しようと試みている状況である。ランダム係数モデルはこのようなケースで用いられることが多い。なお説明変数に連続的な時間 t_{ij}を含めるようなモデルであれば、反復測定を全員が同じタイミングで行っている必要はない。
まずランダム効果を全く考えないことにすると、モデルは以下の単回帰モデル

 \displaystyle
y_{ij} = \beta_0 + \beta_1 t_{ij} + \varepsilon_{ij}

になる。これは、誤差的な変動はあるものの、平均的には全ての個人が最初から最後まで同じ推移をたどることを想定していることになる。

そこでランダム効果モデル、つまり切片だけに個人差(ランダム効果)を考えることにすると

 \displaystyle
y_{ij} = (\beta_0 + r_{i0}) + \beta_1 t_{ij} + \varepsilon_{ij}

という形のモデルになる。ここでは直線の全体的な位置には個人差があることを許容しているが、直線の傾き(測定値の変化するスピード)については全員同じものとみなされている。

さらに、直線の傾き、つまり説明変数である t_{ij}の係数にもランダム効果を設定すると

 \displaystyle
y_{ij} = (\beta_0 + r_{i0}) + (\beta_1 + r_{i1}) t_{ij} + \varepsilon_{ij}

となり、直線の位置も、傾きも各個人で異なるような状況を表現することができる。ランダム係数モデルはこのような形式のモデルとなる。説明変数は、通常の線形モデルと同様に複数含めてもよい。

f:id:mstour:20210220203518j:plain


ところで、前回の記事で述べたように混合モデルは一般に

 \displaystyle
\boldsymbol{y} = X \boldsymbol{\beta} + Z \boldsymbol{u} + \boldsymbol{\varepsilon} \hspace{20pt} (1)

の形式で表すことができる。
上記のランダム係数モデルをどのようにして一般形(1)で書けるかを確認しよう。
 \boldsymbol{y}は目的変数 y_{ij}を全て並べたベクトル、同じく \boldsymbol{\varepsilon}は誤差 \varepsilon_{ij}を並べたベクトルであり、パラメータベクトルは \boldsymbol{\beta} = (\beta_0, \beta_1)^T、ランダム効果のベクトルは \boldsymbol{u} = (r_{10}, r_{11}, r_{20}, r_{21}, \cdots, r_{N0}, r_{N1})^T(全員のランダム効果を全て並べたもの)になる。
またデザイン行列 Xは通常の線形モデルと同様に

 \displaystyle
X = \left(
\begin{array}{cccc}
1 & t_{11} \\
1 & t_{12} \\
\vdots & \vdots \\
1 & t_{1T} \\
1 & t_{21} \\
\vdots & \vdots \\
1 & t_{NT}
\end{array}
\right)

であり、ランダム効果に関するデザイン行列 Z

 \displaystyle
Z = \left(
\begin{array}{cccc}
1 & t_{11} & 0 & 0 & \ldots & 0 & 0 \\
1 & t_{12} & 0 & 0 & \ldots & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
1 & t_{1T} & 0 & 0 & \ldots & 0 & 0 \\
0 & 0 & 1 & t_{21} & \ldots & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
0 & 0 & 0 & 0 & \ldots & 1 & t_{NT} \\
\end{array}
\right)

という形となる。

ランダム係数モデルも式(1)の一般形で表せるので、固定効果パラメータベクトル \boldsymbol{\beta}は最尤法を用いて

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

と推定できる。
また目的変数ベクトル \boldsymbol{y}の分散共分散行列は

 \displaystyle
V(\boldsymbol{y}) = V = Z G Z^T + R

となる。ただし G, Rはそれぞれランダム効果ベクトル \boldsymbol{u}と誤差ベクトル \boldsymbol{\varepsilon}の分散共分散行列である*2

ここまで、一つの説明変数の係数にランダム効果を設定した場合を考えてきたが、全く同じようにして複数の係数にランダム効果を設定することができ、同様に一般形で書き表すことができる。


応用例

今回は、教科書的な説明はそこそこにして、ランダム係数モデルの応用例を紹介したい。
事例として、「全身性強皮症に伴う間質性肺疾患」(SSc-ILD)に対するニンテダニブ(Nintedanib)の治療効果を検討したランダム化比較臨床試験を取り上げる[3][4][5](全身性強皮症(Systemic sclerosis)については例えばこちらなどを参照されたい)。

この臨床試験は、SSc-ILD患者をニンテダニブを投与するグループとプラセボ(偽薬)を投与するグループにランダムに分け、両グループを比較することによってニンテダニブの有効性と安全性を評価したものである。有効性の主要な評価指標は「1年あたりの努力肺活量(Forced Vital Capacity; FVC)の変化量」(以下、参考資料に従い「年間減少率(Annual rate of decline)」*3とする)であり、ニンテダニブを投与することによりFVCの減少が抑えられるかどうかが検討された。
主要な評価のための解析モデルとして、ランダム係数モデルが用いられた。モデルの構造は以下の通りである。なお実際にはもっと細かな設定も必要だが、今の段階では説明不足なので省略する。

[説明変数]
投与群(ニンテダニブorプラセボ)、ATA抗トポイソメラーゼI抗体)の有無、性別、時間、ベースライン時(投与される前)のFVC、年齢、身長、投与群と時間の交互作用、ベースライン時FVCと時間の交互作用
[ランダム効果]
切片、時間の係数

説明変数に時間を含め、ランダム効果を切片および時間に設定しているのは最初に紹介した例と同様であるが、その他にも多数の説明変数が含められている。ここで関心があるのは(横軸を時間としたときの)直線の傾きの治療群ごとの違いである。
時間変数を「年」の単位にすることによって、時間に関する係数(直線の傾き)がFVCの年間減少率の推定値になる(説明変数の係数は、その説明変数が「1単位変化した時の」目的変数の変化量を表していることに注意されたい)。また、「投与群と時間との交互作用」が、ニンテダニブを投与したグループとプラセボを投与したグループとの傾きの差(つまり、FVC年間減少率の差)にあたる。

FVC年間減少率の推定結果は以下の通りであった(資料[3][4]を元に筆者作成)。

f:id:mstour:20210220040310p:plain

この結果の算出手順を整理してみよう。

まず、上記で挙げた説明変数とランダム効果を含むランダム係数モデルをデータに当てはめた。ここで「投与群、時間、投与群と時間の交互作用」以外の説明変数を省略して記載すると、モデルは以下のように表せる:

 \displaystyle
FVC_{ij} = (\beta_0 + r_{i0}) + (\beta_1 + r_{i1}) TIME_{ij} + \beta_2 DRUG_{i} + \beta_{12} TIME_{ij} \times DRUG_i + \varepsilon_{ij} \hspace{20pt} (2)

 iは患者、 jは測定値の通し番号を示すインデックスとする。また「FVC」は反復測定されたFVCの各値、「TIME」は j回目測定の投与開始時点からカウントした時間(単位は年)、「DRUG」は投与群(グループ)でありニンテダニブ投与のグループなら1、プラセボ投与のグループなら0をとるものとする。「TIME  \times DRUG」は投与群と時間の交互作用で、それぞれの変数を掛け算したものである。
このようなランダム係数モデルを用いることによって、FVCの変化には個人差があることを考慮に入れたうえで、集団全体での変化をモデリングしていることになる。

f:id:mstour:20210220203534j:plain

次に、モデルの各固定効果パラメータを最尤法で推定する*4。また、信頼区間の算出や検定を行うために必要となる、パラメータ推定量の分散も推定する(統計ソフトでは通常これらもまとめてアウトプットされる)。

推定結果から、FVCの年間減少率が算出できる。ニンテダニブ群の平均は上記のモデル式(2)のパラメータ \beta_1 \beta_{12}の推定値を足したもの、プラセボ群の平均はパラメータ \beta_1の推定値となる。それぞれの標準誤差は、計算方法が少し複雑なので省略する。この辺りの理論的詳細はまたいずれ書きたい。
2グループの差の平均は、パラメータ \beta_{12}の推定値が該当する*5。差の信頼区間の前後それぞれの幅は \beta_{12}の標準誤差に(通常は)1.96を掛けた値が用いられる。


まとめ

今回は、混合モデルの紹介のところで書ききれなかったランダム係数モデルについて取り上げた。
ランダム係数モデルは、モデルの切片(つまり、目的変数の平均的な水準)だけでなく、説明変数の係数にもランダム効果を設定することによって、説明変数の効果(直線で言うと、傾き)にも個人差があることを考慮したモデルとなっている。
また医学分野での実際の応用例も紹介した*6

今回は以下の資料を参考にした。
[1] Annette J. Dobson(2008), "一般化線形モデル入門", 共立出版.
[2] Brown and Prescott(2015), "Applied Mixed Models in Medicine", Wiley.
[3] Distler et al.(2019), "Nintedanib for Systemic Sclerosis–Associated Interstitial Lung Disease", N Engl J Med.
https://www.nejm.org/doi/full/10.1056/nejmoa1903076
(使用されたランダム係数モデルの詳細はリンク先付録のProtocolを参照。アクセスできない方には申し訳ありませんが、結果だけなら次の資料[4]の9ページ目にも記載されています。)
[4] 医薬品医療機器総合機構, "ニンテダニブエタンスルホン酸塩 審査報告書".
https://www.pmda.go.jp/drugs/2019/P20191223003/530353000_22700AMX00693_A100_1.pdf
[5] Hu et al.(2021), "Power and sample size for random coefficient regression models in randomized experiments with monotone missing data", Biometrical Journal.

*1:一方、Covariance pattern modelのようにランダム効果を導入せずに正規分布以外の確率分布を扱う方法としては、一般化推定方程式(Generalized estimating equation; GEE)がよく知られている。いずれ解説を書いてみたい。

*2:ランダム効果が1個だけのランダム効果モデルとは異なり、ランダム効果が複数あるランダム係数モデルでは、同じ個人のランダム効果どうしの相関を考える必要がある。

*3:このような表現になっているのは、1年あたりの変化という「時間あたりの量」であるため「率(Rate)」であること、またFVCは疾患の進行に伴って減少していくこと、が理由だと思われる。

*4:実際には、制約つき最尤法(Restricted maximum likelihood method; REML)という推定方法を用いることが多い。

*5:端数を丸めているので、上記の表ではニンテダニブ群の平均からプラセボ群の平均を引いたものと差の平均とは少しずれている。

*6:筆者の見聞する限りでは、臨床試験でこのモデルが使われていることはそこまで多くはない。