多施設共同研究では、施設ごとに患者背景や診療内容が異なることがあります。また、同じ施設の患者同士は互いに似た特徴を持つため、たとえ300例の患者がいても、独立した300例分の情報が得られているとは限りません。
「このようなデータの構造を無視すると、標準誤差が過小評価され、信頼区間が実際より狭くなってしまう」よく説明されますが、ただし、これは必ずしもいつでも起こるわけではありません。影響の大きさは、関心のある変数が「クラスター（施設）の中でどの程度ばらつくのか」、あるいは「施設単位で決まるのか」によって変わります。
効く場合と効かない場合
2つの状況をシミュレーションで比べます。40クラスター、各10例、クラスター効果の標準偏差0.8。片方は曝露が個人単位で決まる状況、もう片方はクラスター単位で決まる状況です。以下は説明のための架空データです。
曝露が個人ごとに決まるなら、比較はクラスターの内側で行われます。施設の予後の良し悪しは治療群と対照群の両方に等しくかかるので、群間差の推定にはほとんど影響しません。クラスタリングの調整は要らないとまでは言いませんが、効果は小さくなります。
一方、曝露がクラスター単位で決まる場合は話が別です。たとえば、クラスター無作為化試験、施設単位の診療方針、地域単位の介入などがこれに当たります。
この場合、実質的な標本数を決めるのは患者数ではなく、クラスターの数です。たとえば40クラスター、400例のデータを、クラスター構造を無視して「400例の独立したデータ」として扱えば、標準誤差を過小評価し、信頼区間を不当に狭くしてしまいます。極端に言えば、40クラスターを400例として扱うことで、信頼区間の幅が4割程度まで狭くなることもあります。
そのため、データを受け取ったときに最初に確認するのは、関心のある変数がクラスターの内側で変動しているかどうかです。ここが決まれば、対処の必要性も決まります。
ハザードに個体差を掛ける
frailtyモデルは、クラスターごとの予後の違いをハザードの乗数として表します。
$$h_{ij}(t) = z_i \, h_0(t)\exp(\boldsymbol{\beta}^\top \mathbf{x}_{ij}) \tag{1}$$
式(1)の$z_i$がクラスター$i$のfrailtyです。$z_i$が2なら、そのクラスターの患者は基準の2倍のハザードを持ちます。観測されない施設要因をまとめて1つの乱数に押し込んだ形で、線形混合モデルの変量切片に相当します。
$z_i$には平均1の分布を置きます。ガンマ分布が標準で、平均を1に固定すると分散1つで形が決まります。
この分散が、frailtyモデルから得られる固有の情報です。ロバスト分散は標準誤差を直すだけで、クラスター間のばらつきの大きさは教えてくれません。施設間差そのものが研究の関心事なら、frailtyモデルを選ぶ理由になります。
survivalパッケージのrats（100腹×3匹、42イベント）で4通りを比べます。
coxph(Surv(time, status) ~ rx, data = rats) # 無視
coxph(Surv(time, status) ~ rx + cluster(litter), data = rats) # ロバスト分散
coxph(Surv(time, status) ~ rx + frailty(litter), data = rats) # frailty
coxph(Surv(time, status) ~ rx + strata(litter), data = rats) # 層別
推定されたfrailtyの分散は2.02でした。腹による差はかなり大きいことになりますが、それでも治療効果の推定値と区間は4通りでほとんど変わりません。クラスター間差が大きいことと、治療効果の推定が歪むことは別の話です。
3つの選択肢の使い分け
ロバスト分散は、点推定をそのままに標準誤差だけを直します。モデルの形を変えないので、報告されるハザード比の意味は通常のCox回帰と同じ、集団平均的な効果です。クラスター数が十分にあれば（目安として40程度以上）、これがいちばん手数が少なく済みます。
frailtyモデルは、クラスター間のばらつきを推定します。得られるハザード比は「同じ施設の中で比べたときの効果」という条件付きの意味になり、ロバスト分散のものとは解釈が異なります。frailtyの分散が大きいほど両者は離れます。施設間差の大きさを報告したい場合や、クラスター数が少なくロバスト分散が不安定な場合に選びます。
層別Cox回帰では、クラスターごとに基準ハザードを別々に設定します。基準ハザードについて強い分布仮定を置かなくてよいことが利点ですが、その分、クラスターごとに情報を分けて使うため、クラスターが小さいと情報を大きく失います。1施設あたり数例しかないような研究では、あまり勧められません。
私が多施設研究の相談を受けたときは、まず施設数と1施設あたりの症例数を聞きます。施設数が多くて症例が散らばっているならロバスト分散、施設数が少ないか施設間差そのものを論じたいならfrailty、施設が大きくて数が限られているなら層別も候補になります。
クラスター数が少ないとき
ロバスト分散は大標本の性質に頼っています。ここでの大標本はクラスター数で、症例数ではありません。クラスターが10や15しかない研究では、ロバスト分散が真の分散を過小に推定することが知られています。
目安として、クラスター数が30を下回るなら注意が要ります。対処としては、自由度をクラスター数から調整したt分布で区間を作る方法や、クラスター単位のブートストラップがあります。frailtyモデルを使う選択もありますが、こちらはクラスター数が少ないと分散パラメータ自体が不安定になるので、万能ではありません。
もうひとつ、frailty分布の指定にも仮定が入ります。ガンマ分布が既定で使われるのは計算が扱いやすいためで、生物学的な根拠があるわけではありません。対数正規分布を指定することもできます。分散の推定値は分布によって変わるので、この値を主要な結果として報告するなら、他の分布でも当てて感度を見るのが安全でしょう。
施設ごとの予後を推定する
frailtyモデルは、施設ごとの$z_i$の事後的な推定値も返します。これを並べれば、どの施設の予後が良いかが分かります。医療の質の評価に使われる形です。
ここで効くのが縮小です。症例数の少ない施設では、推定値が全体平均のほうへ引き寄せられます。3例しかない施設でたまたま全員が亡くなっても、その施設のfrailtyが極端に大きく推定されることはありません。単純に施設ごとの死亡率を並べるより、順位が安定します。
この性質は、施設間比較を公表する場面で重要になります。縮小のない生の率で順位をつけると、上位も下位も小規模施設が占めることが知られています。偶然の振れ幅が大きいためです。
ただし、縮小には注意すべき点もあります。本当に成績の悪い小規模施設であっても、推定値は全体平均のほうへ引き寄せられるため、問題の検出が遅れる可能性があります。
したがって、施設の順位だけを示すよりも、各施設の推定値とその信頼区間を示し、全体平均から明確に区別できる施設なのかどうかを見るほうが安全です。施設間比較では、「順位をつけること」と「統計的に明らかな差があること」を分けて考える必要があります。
frailtyと時変効果を混同しない
frailtyが存在すると、集団全体で観察したハザード比が時間とともに1へ近づいていくことがあります。これは、frailtyの大きい、つまり脆弱な個体ほど早くイベントを経験して集団から脱落し、時間の経過とともに残っている集団の構成が変化するためです。
このような見かけ上の時間変化は、治療効果そのものが時間とともに減弱している場合と区別しにくい点に注意が必要です。たとえば、Schoenfeld残差を用いた検定で比例ハザード性の違反が認められたとしても、それだけでは、その原因が治療効果の減衰なのか、未観測の異質性によるものなのかをデータから判別することはできません。
そして、両者では臨床的な意味が大きく異なります。治療効果が時間とともに減衰しているのであれば、「治療効果は長期的には弱くなる」と解釈します。一方、異質性によって見かけ上ハザード比が1へ近づいているのであれば、「個々の患者における治療効果は持続しているが、脆弱な患者が先にイベントを経験することで、集団全体で観察した効果が時間とともに弱く見えている」と考えられます。
frailtyモデルでクラスター間・個体間の異質性が大きいことが示された場合には、観察された時間依存性について、真の治療効果の変化だけでなく、未観測の異質性による選択の影響も考慮する必要がある、という考察につなげることができます。
報告の仕方
クラスターの単位と数、1クラスターあたりの症例数の分布を書きます。「12施設、1施設あたり8-64例（中央値23例）」のような記述です。この情報がないと、読者は分散推定の妥当性を判断できません。
どの方法で対処したかも明記します。ロバスト分散ならその旨、frailtyなら分布と推定された分散、層別なら層別因子。frailtyの分散は載せる価値があります。0に近ければクラスタリングは実質的に問題なかったことになり、大きければ施設間差が大きいという知見になります。
そして、関心のある変数がクラスター内で変動しているかどうかを1文で触れます。図1のとおり、ここが結論の頑健性を左右します。クラスター単位の曝露なら、クラスター数が実質的な標本数であることも書いておくと、査読者との議論が短くて済みます。
frailtyの分布をどう選ぶか
ガンマ分布が既定で使われるのは、周辺尤度が閉じた形で書けて計算が軽いためです。
また、対数正規分布も選べます。線形混合モデルの変量効果と同じ形なので、他の解析との接続が良好になりますが、推定は数値積分が入るぶん重くなります。
そして、分散の推定値は分布によって変わります。ratsの例ではガンマで2.02でしたが、対数正規を当てれば別の値になります。分散そのものを主要な結果として報告するなら、両方で当てて値の幅を並べるほうが誠実でしょう。
治療効果の推定値は、分布の選択にそれほど敏感ではありません。図3で4通りの扱いが2.0前後に収まったのと同じで、関心のある係数への影響は限定的です。分布の選択が効くのは、クラスター間差そのものを論じるときです。
次に確かめること
反復イベントのデータも同じ構造を持ちます。1人の患者に複数回のイベントがあれば、患者がクラスターです。この場合はfrailtyモデルかロバスト分散が要ります。
層別との比較は層別Coxモデル、モデルの基本形はCox比例ハザードモデル、比例ハザード性の確認はSchoenfeld残差にあります。frailtyを入れても比例ハザード性の仮定は残るので、確認は必要です。
施設をどう扱うか、クラスター数が少ないときに何を選ぶか。このあたりは研究の規模と施設あたりの症例数を見ないと決まりません。Dr.データサイエンスでは、こうした解析設計のご相談を承っています。
参考文献
Hougaard P. Frailty models for survival data. Lifetime Data Analysis. 1995;1(3):255-273.
Lin DY, Wei LJ. The robust inference for the Cox proportional hazards model. Journal of the American Statistical Association. 1989;84(408):1074-1078.
Glidden DV, Vittinghoff E. Modelling clustered survival data from multicentre clinical trials. Statistics in Medicine. 2004;23(3):369-388.
Duchateau L, Janssen P. The Frailty Model. Springer; 2008.
Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000.
