« 読了：Navarro, et al. (2005) ディリクレ過程に基づく個人差モデリング | メイン | 読了：Scott, McPherson, Ramsay, Campbell (2002) ランダム化比較試験でたくさんの共変量を条件間でバランスさせる割付方法 »
2016年8月 3日 (水)
Zubizarreta, J.R. (2012) Using Mixed Integer Programming for Matching in an Observational Study of Kidney Failure After Surgery. Journal of the American Statistical Association, 107(500), 1360-1371.
ケース・コントロール・スタディで共変量が多いとき、対象者のマッチングを混合整数計画(MIP)でやりましょうという論文。ネットに落ちてたdraftで読んだ。
まず、観察研究におけるマッチングとはなにか、その利点について説明があって...
マッチングによくつかわれているのはネットワーク・アルゴリズムである。これは処理単位とそれにマッチングするコントロールのあいだの距離の合計を反復で最小化させる方法だが、どれだけ反復すればいいのかわからないし、コルモゴロフ・スミルノフ検定料とか共変量間の相関といった統計量を明示的にバランシングすることができない。
そこで新手法をご提案しますよ。Rのパッケージmipmatchもつくりましたよ。
この手法ではMIPを使います。なぜこれまで使われてこなかったか。それはtractableでなかったから(所与のサイズのもっとも難しい問題を解くのに必要なステップ数がサイズの多項式になっていなかったから)。これに対してネットワーク・アルゴリズムなら、シンプレクス法とかオークション・アルゴリズムとかで解けた。Rのoptmatchを参照のこと。
しかし[...と線形計画についての紹介があって...]、最近MIPはとても発展し、マッチングに使えるようになったのであります。
処理単位の集合を$T=\{t_1, \ldots, t_T\}$、潜在的なコントロール単位の集合を$C=\{c_1, \ldots, c_C\}$としよう($T \leq C$)。共変量のラベルの集合を$P = \{p_1, \ldots, p_P\}$とする。任意の処理単位とコントロール単位について、それぞれが共変量のベクトル$x_t, x_c$を持ち、距離$0 \leq \delta_{t,c} < \infty$を定義できるとしよう。
さあ、整数計画として定式化しまっせ。深呼吸...
変数$a_{t,c}$を、処理単位$t$がコントロール単位$c$に割り当てられたときに1、どうでないときに0とする。ベクトル(というか行列)にまとめて$a$とする。
以下の式を最小化する$a$を
$\sum_{t \in T} \sum_{c \in C} \delta_{t,c} a_{t,c} + \sigma_{i \in I} w_i \mu_i (a)$
第一項はマッチドペアの共変量の総和。第二項では、共変量のインバランスを表すいくつかの指標$mu_i (a)$があり、それに重み$w_i$をつけて足しあげている。指標をどうするかはあとで述べる。
制約条件は次の4本。
1本目、$\sum_{c \in C} a_{t, c} = m, t \in T$。それぞれの処理単位は $m$ 個のコントロール単位とマッチする。
2本目、$\sum_{t \in T} a_{t, c} \leq 1, c \in C$。それぞれのコントロール単位は、0個ないし1個の処理単位とマッチする。
3本目、$a_{t,c} \in \{0,1\}, t \in T, c \in C$。そりゃそうですわね。
4本目、$v_j(a) \leq e_j, j \in J$。ここで$v_j(a)$とは、制約式の$\mu_i(a)$と同様に、共変量のインバランスを表すなんらかの指標。[←なるほどね。制約をかける指標と最適化する指標が違っていていいように別の記号を使ってるわけだ]
では、$\mu_i(\cdot)$をどう定義するか。具体的事例に即して考えます。[以下、論文の表記よりも簡略化してメモする]
その１、単変量のモーメントをバランスさせたいとき。
$\mu_i(a) = | \sum_t \sum_c \frac{x_{c,i} a_{t,c}}{mT} - \bar{x}_{T,i}|$
ここで$i$は共変量のインデクス。$\bar{x}_{T,i}$とはすなわち、処理群における共変量$i$の平均。$\sum_c \frac{x_{c,i} a_{t,c}}{m}$で、処理単位$t$にマッチした$m$個のコントロール単位における共変量$i$の平均になるから、つまり処理群における平均からの各コントロール単位の偏差の絶対値の和を最小化しているわけだ。[←あー、なるほどね。ここでは処理単位のセットが事前に決まっているわけね]
絶対値記号が入ってて線形になっていないじゃんとお思いのみなさん。こうするといいのだよ。[←線形計画定式化の必須テクである。文系の私には、これが腑に落ちるまで結構な時間がかかった...]
最適化式は
$\sum_{t \in T} \sum_{c \in C} \delta_{t,c} a_{t,c} + \sigma_{i \in I} w_i z_i$
ってことにしといて、制約式に次の2本を追加する。
$z_i \geq \mu_i(a)$
$z_i \geq -\mu_i(a)$
実データを使って実験してみる。$w_i$を大きくすれば満足いく結果が得られる。一般的アドバイスとしては、距離行列を平均で割っておき、共変量も標準化した状態で、処理単位数の1/10の値にするとよい。
このように最適化式でバランスさせるのと制約式でバランスさせるのとどちらがよいか。制約式はわかりやすいし、MIPじゃなくて純粋な整数計画になるという利点がある。いっぽう閾値を小さくし過ぎるとinfeasibleになるというのが欠点。
その２、共変量の多変量のモーメント(相関)をバランスさせる場合には...[略]
その３、コルモゴロフ-スミルノフ検定料をバランスさせる場合には...[略]
その４、共変量を正確にマッチングさせたい場合。
$x_{.,p}$を名義共変量とし、とる値を$b \in B$とする。正確なマッチングは以下の制約式で表現できる。すべての$b$について
$\sum_t \sum_c a_{t,c} I(x_{t,p}=b AND x_{c,p}=b) = m \sum_t I(x_{t,p} = b)$
これはぴったり同じ場合だが、差の絶対値の上限を決めて緩く制約しても良い。なお、当該共変量で分割して各部分でマッチングさせればそりゃぴったり合うけれど、そのやりかただと緩く制約することができない。
とこで、上記のようなマッチングをfine balanceという。いっぽう
$\sum_t \sum_c a_{t,c} I(x_{c,p}=b) = m \sum_t I(x_{t,p} = b)$
というのをnear-fine balanceという。
ネットワーク・アルゴリズムでは複数の名義共変量の周辺分布をfine balanceないしnear-fine balanceさせることはできないが[←その理由が説明してあるが理解できなかった。Rosenbaum(2010)を読めとのこと]、提案手法だとできる。
というわけで、以上を実現するRパッケージ mipmatch を作りました。これは中でIBMのCPLEXを呼びます。[...説明略]
実データへの適用例 ... [略]
$m$が処理単位のあいだで変動してよい場合には ... [略]
考察。
いまやMIPはすごい。使いましょう。
マッチングは傾向スコアでやるという手もある。しかしそれは確率的なバランスが達成できるだけだし、ケース・コントロール研究においては傾向スコアが存在しない[←え、なんで？ なくはないんじゃないの？]
云々、云々。
。。。いやー、この論文も勉強になりましたです。
しっかし、この定式化だと、$a$は(処理単位数 x 潜在コントロール数)の二値行列になり、各行の和が$m$となるという制約と各列の和が0か1になるという制約しかかからないわけだ。サンプルサイズが大きいと、いかにCPLEXとはいえ、手に負えなくなるだろうな...
この話は観察研究におけるマッチングの話だけど、ランダム化比較試験でもいくつかの共変量によって層別することはあるわけで(stratified randomization)、共変量がすごく多ければこの論文と同じ問題を抱えることになる。そういうときにMIPを使っている事例を探しているのだけれど、見当たらない。たぶん医学統計じゃなくて工業試験みたいな分野を探せばいいんだろうけど、検索のキーワードがわかんないんだよなあ。
論文：データ解析(2015-) - 読了：Zubizarreta(2012) ケース・コントロール研究でのマッチングを混合整数計画でやりましょう