第 5 章 Selection Model 40
5.3 MAR を仮定した SM
同じSMであっても,欠測メカニズムとしてMARを仮定するかMNARを仮定するかによって,解析の理 論・実装の容易さや結果の評価の仕方が大きく異なる.以下,この違いにも注目しつつ,式を追いながら解析 手法を検討していく.なお,パラメータが分離していることを仮定しているため,MARの場合は2章で定義 した「無視できる欠測メカニズム」となる.このとき,計算が極めて容易になることを以下で確認する.
5.3.1 MAR を仮定した SM その 1 ( SM の一般論)
MARを仮定し,興味あるパラメータが応答変数の分布のパラメータθのみとする.このとき,Riの(条件 付き)分布は
f(Ri|Xi,Yi,bi,ψ) =f(Ri|Xi,Yio,ψ)
3記号が煩雑となるため,繰り返しておくが
θ= ( θ1
θ2
) である.
となる.これより,完全データの尤度(変量効果あり)は
f(Yi,Ri,bi|Xi,θ,ψ) =f(Yi|Xi,bi,θ1)·f(Ri|Xi,Yio,ψ)·f(bi|Xi,θ2) (5.2) となる.しがたって,観測データの尤度は
f(Yio,Ri|Xi,θ,ψ) =
∫ ∫
f(Yi,Ri,bi|Xi,θ,ψ)dbidYim
=
∫ ∫
f(Yi|Xi,bi,θ1)·f(Ri|Xi,Yoi,ψ)·f(bi|Xi,θ2)dbidYim (∵(5.2))
=f(Ri|Xi,Yio,ψ)·
∫ ∫
f(Yi|Xi,bi,θ1)·f(bi|Xi,θ2)dbidYim
=f(Ri|Xi,Yio,ψ)
∫
f(Yio,Ymi |Xi,θ1,θ2)dYmi
=f(Ri|Xi,Yio,ψ)·f(Yio|Xi,θ) となる.ここで,両辺の対数をとると,
log{f(Yio,Ri|Xi,θ,ψ)}= log{f(Ri|Xi,Yoi,ψ)·f(Yio|Xi,θ)}
= log{f(Ri|Xi,Yoi,ψ)}+ log{f(Yio|Xi,θ)}
となる.つまり,応答変数の分布のパラメータθのみが興味の対象であり,推定に最尤法を用いる場合,脱落 確率の(条件付き)分布f(Ri|Xi,Yoi,ψ)がどのように変わっても,θで微分すると0となるため,θの推定 には影響を与えない.つまり,欠測メカニズムが「無視できる」わけである.
【Yiの分布の仮定】
いま,応答変数Yiの分布については特に規定することなく,一般論で検討してきた.しかし,実際の解析 では,この分布に対して何らかの仮定をすることが普通である.
現実的によく用いられる仮定の1つは,Yiに線形混合モデルを,変量効果・誤差共に正規分布を仮定する ものである.このとき,解析はたとえばSASのPROC MIXED等を用いて容易に行える.一方,正規分布を当 てはめることが妥当であったか,という点に検討の余地が残る.
別の方法として,Yiの分布に極力仮定をしない,という方法も考えられる.というのも,SMでは,f(Yi|Xi,bi) = f(Yio,Ymi |Xi,bi)に分布を仮定しても,欠測がある場合は一般にその分布の仮定が正しいかどうかを確認す ることができない4ので,「できるだけ分布の仮定を少なくすることで,誤特定のリスクを減らしたい」と考え るためである.このような目的で用いられる方法に,IPW,GEE,wGEEやDoubly Robustな推定量等を用い るセミパラメトリックな方法がある5.これらについては8章で述べるが,モデルの仮定が少ないため,誤特定 が減らせる反面,推定量の分散の計算やソフトウェアでの実行がやや面倒となる.また,パラメトリックなモ デルが正しい場合には,推定量の分散をより大きく推定してしまうことが多くなる.
5.3.2 MAR を仮定した SM その 2 ( MMRM )
次に,MARを仮定した一般のSMとの違いに注目しながら,MMRM (Mixed Model for Repeated Measures) について少し触れる.MMRMは9章でより詳しく扱う.
"MMRM"という用語は,文献によって異なる定義で用いられており,モデルや推定方法が厳密に表現されて
いるとは言えない.Mallinckrodt (2013)の7.2節には,「"MMRM"という用語自体の使用が推奨されておらず,
Direct Likelihood (DL)法と呼んだ上で,解析計画を立てる際はモデルの詳細を追記するのが望ましい」旨が述
べられている.一方,本報告書では「モデル・推定方法に対して適切な指定をした上であれば,用語を使用す
4感度分析の部分で詳しく述べるが,Type (i)の検証不能な仮定である.
5NRC (2010)のp73では,このような理由から,5章の感度分析でSMに対してセミパラメトリックなモデルを用いて検討する旨が記
載されている.
ること自体は許容する」という立場をとる.本章ではMMRMという名前のもとに比較的よく用いられている モデル・推定方法に対して,一般のSMの中で特徴的な部分について述べる.
まず,固定効果は「投与群」「時点(カテゴリ値)」「時点と投与群の交互作用」とする.次に,分散共分散 構造を考えるが,線形混合モデルでは変量効果の分布と誤差の分布を解析者が特定する必要がある(詳しくは
Appendix A参照).この部分の指定が,MARを仮定した一般のSMの中で,MMRMの特徴的な部分である.
【MARを仮定した一般のSM】
MARを仮定した一般のSMでは,上でも述べた通り一般に応答変数の分布を特定していない6.しかし,以 下,特にMMRMとの対応を考えるため,応答変数に(多変量)正規分布を仮定したモデルを扱う.
このとき,変量効果と誤差のそれぞれに対して分散共分散行列を指定する.たとえば,3時点,変量効果は 被験者のみ切片項に含め,誤差は独立同分布で分散σ2のとき,
Yi=Xiβ+bi13+ϵi (biiid∼ N(0, ν2), ϵiiid∼N(0, σ2I3), biとϵiは独立) となる.このとき,周辺分散は
V[Yi] =V[bi13+ϵi]
=13V[bi]1′3+V[ϵi]
=13ν21′3+σ2I3
=ν2
1 1 1 1 1 1 1 1 1
+σ2
1 0 0 0 1 0 0 0 1
=
ν2+σ2 ν2 ν2 ν2 ν2+σ2 ν2 ν2 ν2 ν2+σ2
となる.なお,変量効果として被験者を切片項に入れる場合,通常は誤差分散に無構造(unstructured)を用いる ことはない.理由は以下のMMRMの項で説明する.
【MMRM】
MMRMでよく用いられるモデルでは誤差分散を無構造(unstructured)とする.ここでたとえば,3時点,変 量効果は被験者のみ切片項に含め,誤差分散をunstructuredとし,変量効果と誤差が独立とすると
Yi=Xiβ+bi13+ϵi
bi
iid∼N(0, ν2), ϵi iid∼N
0 0 0
,
σ21 σ12 σ13
σ12 σ22 σ23 σ13 σ23 σ23
となる.このとき,周辺分散は
V[Yi] =V[bi13+ϵi] (5.3)
=13V[bi]1′3+V[ϵi] (5.4)
=13ν21′3+
σ12 σ12 σ13 σ12 σ22 σ23
σ13 σ23 σ23
6本報告書の対象外であるため深入りしないが,カテゴリカルデータに対して一般化線形混合モデルを用いることも可能である.
=ν2
1 1 1 1 1 1 1 1 1
+
σ21 σ12 σ13
σ12 σ22 σ23
σ13 σ23 σ23
=
ν2+σ21 ν2+σ12 ν2+σ13
ν2+σ12 ν2+σ22 ν2+σ23 ν2+σ13 ν2+σ23 ν2+σ23
(5.5)
となり,全ての成分にν2が出てくる.いま,異なる被験者のデータは全て独立のため,この状況ではν2, σ21, σ12,· · · の推定値を1つに絞ることができない7.
以上より,変量効果と誤差(相関構造:unstructured)の両方を指定した場合,分散成分のパラメータを1つ 1つ推定することはできないことがわかる8.また,「変量効果を特定する」ということは,「誤差の持っている 複雑でよくわからない時点間相関・不等分散などを,変量効果という解釈可能な成分に特定(モデリング)し てやり,その分だけ誤差の分布をより簡単なものにする」こととも考えることができる.そのように考えるな らば,変量効果を指定した後でも誤差分散がunstructured,ということは,変量効果を指定しても誤差の複雑さ が改善されなかった,ということとなる.
MMRMという手法が提案されたMallinckrodt et al. (2001)では,MMRMのモデルの式中に変量効果として 被験者を用いている.その上で,誤差分散をunstructuredとして,推定の際は変量効果は誤差と合わせて推定 している.これは一見奇妙なようにも思えるが,以下のように解釈できる.
• 1被験者の経時データは全て相関していると考えられるため,変量効果として被験者を考慮することは 自然である.
• その上で,誤差の分散共分散構造は事前に想定できないため,unstructuredを指定する.
• このとき,解析する際に,変量効果と誤差を別々に指定すると,パラメータの推定ができなくなるた め,「変量効果の分散+誤差分散」 を1パラメータとして考える.具体的には,5.5)中の”ν2+σ12”や
”ν2+σ12”のような「足し算されたもの」をそれぞれ1パラメータと考え,推定する.
つまり,モデルに変量効果は指定するものの,全てのパラメータの推定はできないため,変量効果の分散と合 わせたものだけを推定する,と考えるのである.
では,次にどうしてこのような面倒なことを考えるのかについて,Mallinckrodt (2013) 7.4節に即して述べる.
まず,事前情報の量にもよるが,相関構造を事前に予測することは難しいことが多いこと,実際にunstructured が最もデータに適合することが多いことが挙げられている.さらに,Mallinckrodt et al. (2008)での検討による と,変量効果を正しく特定するのではなくこのように扱ったとしても,薬剤の効果の推定にはほぼ影響を与え ないようである.また,Mallinckrodt et al. (2004)での検討によると,最尤法を用いて推定した結果,相関構造 の誤特定によりαエラーのインフレが起こることがあるようである.これらの観点から,MMRMでは上記の ようなモデリングが基本となっているようである.一方で,血圧や臨床検査値等のようにデータの被験者内相 関のうち,被験者特有の要因が最も大きいような場合には,Compound Symmetryのような構造が適切である かもしれない.このような倹約な(parsimonious)モデルは,(正しく特定できるならば)unstructuredを仮定した モデルより検出力が上がる,等のメリットがある.常にunstructuredの使用が推奨されているわけではない点 にも,注意が必要である. また,相関構造を1つに特定するのではなく,複数の候補を設定し,情報量規準等 でモデル選択を行うと同時に推定・検定する方法なども考えられる.その場合,モデル選択と推定・検定を同
7識別性(identifiablity)がない,という.具体例で考えると,
ν2+σ21 ν2+σ12 ν2+σ13
ν2+σ12 ν2+σ22 ν2+σ23
ν2+σ13 ν2+σ23 ν2+σ32
=
1 0.9 0.8 0.9 1 0.9 0.8 0.9 1
の場合,例えば
1. ν2= 0.8でσ12=σ22=σ23= 0.2, σ12=σ23= 0.1, σ13= 0 2. ν2= 0.5でσ12=σ22=σ23= 0.5, σ12=σ23= 0.4, σ13= 0.3
の区別をつけることができない.0< ν2<1の他の場合も同様に区別がつかない.
8ただし,変量効果として施設等,被験者以外の要因を追加することは可能である.
一データで行うことの妥当性について,判断が必要となるであろう.妥当でないと判断された場合には,相関 構造に順番をつけておき,最初に収束した構造を選択する,という方法もありうるかもしれない.これらの点 については,その正当性や統計的性能について,引き続き検討が必要であろう.
5.3.3 MNAR (outcome-dependent dropout) を仮定した SM
次に,欠測のメカニズムとしてMNARのうち,Little (2008)で定義されたoutcome-dependent dropoutを仮定 した場合を考える9.つまり,
f(Ri|Xi,Yi,bi,ψ) =f(Ri|Xi,Yi,ψ)
を仮定する.
このとき,完全データの尤度は
f(Yi,Ri,bi|Xi,θ,ψ) =f(Yi|Xi,bi,θ1)·f(Ri|Xi,Yi,ψ)·f(bi|Xi,θ2) となる.ここで,まず変量効果biで積分すると
f(Yi,Ri|Xi,θ,ψ) =
∫
f(Yi|Xi,bi,θ1)·f(Ri|Xi,Yi,ψ)·f(bi|Xi,θ2)dbi
=f(Yi|Xi,θ)·f(Ri|Xi,Yi,ψ) (5.6) となる.(5.6)で変量効果biが脱落確率モデル(f(Ri|Xi,Yi,ψ))に含まれていないことで,積分の計算が可 能となっている10.これを用いて,観測データの(周辺)尤度は
f(Yio,Ri|Xi,θ,ψ) =
∫
f(Yoi,Yim|Xi,θ)·f(Ri|Xi,Yoi,Ymi ,ψ)dYim (5.7) と表せる.MARの場合と異なり,一般に積分はこれ以上計算できない.そこで,Diggle and Kenward (1994)に ならい,欠測識別変数Riに以下のモデルを仮定して,積分の数値計算を行う.Rij = 1が「観測」,Rij = 0 が「欠測」であることに注意して,脱落確率モデルを
logit{P r(Rij = 0|Ri1= 1,· · · , Ri,j−1= 1,Yi,Xi,ψ)}=ψ1+ψ2Yi,j−1+ψ3Yij (j= 2,· · ·, n) (5.8) P r(Rij= 0|Ri,j−1= 0,Yi,Xi,ψ) = 1 (j= 2,· · ·, n) (5.9)
P r(Ri1= 1|Yi,Xi,ψ) = 1 (5.10)
とおく11.このとき,ψ= (ψ1, ψ2, ψ3)′である.また,式(5.8)は,「被験者iの時点jのデータYijの脱落に影 響を与えているのは,1時点前の値Yi,j−1とその時点の値Yijの高々2つである.特にψ3= 0のとき,脱落は 観測値Yi,j−1のみに影響を受けるため,MARとなる」ことを意味し,式(5.9)は「単調な欠測」を意味する.
式(5.10)は,「ベースライン(j= 1)の値は必ず測定される」ことを意味する.
なお,このモデルは以下のことなどを仮定している.
9再度注意しておくが,本節の内容は一般には応答変数に正規性を仮定していない.たとえば,NRC (2010)の5章の感度分析では,応 答変数に正規性は仮定せず,セミパラメトリックなモデルに対してIPW法を用いて推定する方法が提示されている.
10ここで,f(Yi|Xi,bi,θ1)とf(Yi|Xi,θ)の違いは,たとえば線形混合モデル Yi=Xiβ+Zibi+ϵi, bi
iid∼ N(0,D),ϵi
iid∼N(0,Σ)
で{bi}と{ϵi}も独立の場合,
f(Yi|Xi,bi,θ1) =⇒ Yi|bi∼N(Xiβ+Zibi,Σ),
f(Yi|Xi,θ) =⇒ Yi∼N(Xiβ,ZiDZ′i+Σ)
という違いである.ここで,θ2は変量効果biの分散共分散行列Dの成分のパラメータベクトルであり,Yi|biの分布には含まれてい ないが,Yiの分布には含まれていることにも注意が必要である.
11左辺に出てくるXiが右辺では消えているのは,Xiに依存しないモデルであることを意味している.