比例ハザードモデル

比例ハザードモデル（ひれいハザードモデル、英: proportional hazards models）は、統計学における生存モデルの一種である。生存モデルは、ある事象が発生する前の経過時間を、その時間量に関連する可能性がある1つまたは複数の共変量（英語版）に関連づけるものである。比例ハザードモデルでは、共変量の単位増加による固有の効果は、ハザード率に対して乗法的に作用する。たとえば、ある薬を服用すると、脳卒中発症のハザード率が半分になったり、製造部品の構成材料を変更すると、故障のハザード率が2倍になったりする。加速故障時間モデル（英語版）などの他の種類の生存モデルは、比例ハザードを示さない。加速故障時間モデルは、ある事象の生物学的または機械的な生活史が加速（または減速）される状況を記述する。 生存モデルは、2つの部分から構成されていると見なすことができる。基本となるベースライン・ハザード関数（ λ 0 ( t ) {\displaystyle \lambda _{0}(t)} と表記されることが多い）は、ベースラインレベルの共変量で、時間単位あたりのイベントリスクが時間とともにどのように変化するかを記述するものである。一方、効果パラメータは、説明共変量に応じてそのハザードがどのように変化するかを記述するものである。典型的な医療事例では、交絡の変動を抑制し（または）制御するために、治療割り当てや、試験開始時の年齢、性別、および試験開始時の他の疾患の有無など、患者特徴の共変量が含まれる。 「比例ハザード条件」（proportional hazards condition）とは、共変量がハザードに乗法的に関係していることを意味する。定常係数の最も単純なケースでは、たとえば、ある薬剤による治療により、任意の時間 t {\displaystyle t} における被験者のハザードが半減し、一方でベースラインのハザードは変化する可能性がある。ただし、これによって被験者の寿命が2倍になるわけではないことに注意すること。寿命に対する共変量の正確な効果は、 λ 0 ( t ) {\displaystyle \lambda _{0}(t)} の種類に依存する。共変量はバイナリ予測変数に限定されるものではなく、連続的な共変量 x {\displaystyle x} の場合、一般的にハザードは指数関数的に応答すると仮定される。 x {\displaystyle x} の単位増加ごとにハザードは比例的にスケーリングされる。 以下に示すCoxの部分尤度は、ベースライン・ハザード関数のBreslow推定量を用い、それを完全な尤度に当てはめて、その結果が2つの要因の積であることを観測することで得られる。第1の因子は、ベースライン・ハザードが「相殺」された後述の部分的な尤度である。第2の因子は、回帰係数を含まず、打ち切りパターンを介してのみデータに依存する。このように、比例ハザードモデルによって推定された共変量の効果は、ハザード比として知ることができる。 デイヴィッド・コックス卿（英語版）は、比例ハザードの仮定が成り立つ（あるいは成り立つと仮定される）場合、ハザード関数を考慮せずに効果パラメータを推定することが可能であると述べた。このような生存データへのアプローチは、Cox比例ハザードモデルの適用と呼ばれ、Coxモデルまたは比例ハザードモデルと略されることもある。ただし、Coxは、比例ハザード仮定（proportional hazards assumption）の生物学的解釈が非常に難しい場合があることも言及した。 被験者 i の共変量の実現値を Xi = (Xi1, … , Xip) とする。Cox比例ハザードモデルのハザード関数は、 λ ( t | X i ) = λ 0 ( t ) exp ⁡ ( β 1 X i 1 + ⋯ + β p X i p ) = λ 0 ( t ) exp ⁡ ( X i ⋅ β ) {\displaystyle \lambda (t|X_{i})=\lambda _{0}(t)\exp(\beta _{1}X_{i1}+\cdots +\beta _{p}X_{ip})=\lambda _{0}(t)\exp(X_{i}\cdot \beta )} の形式である。この式は、共変量ベクトル（説明変数）Xi を持つ被験者 i の時刻 t におけるハザード関数を示している。 時刻 Yi において、観測されるべき事象が被験者 i に発生する可能性は、次のように書くことができる。 L i ( β ) = λ ( Y i ∣ X i ) ∑ j : Y j ≥ Y i λ ( Y i ∣ X j ) = λ 0 ( Y i ) θ i ∑ j : Y j ≥ Y i λ 0 ( Y i ) θ j = θ i ∑ j : Y j ≥ Y i θ j {\displaystyle L_{i}(\beta )={\frac {\lambda (Y_{i}\mid X_{i})}{\sum _{j:Y_{j}\geq Y_{i}}\lambda (Y_{i}\mid X_{j})}}={\frac {\lambda _{0}(Y_{i})\theta _{i}}{\sum _{j:Y_{j}\geq Y_{i}}\lambda _{0}(Y_{i})\theta _{j}}}={\frac {\theta _{i}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}} ここに、θj = exp(Xj ⋅ β) であり、その合計は時刻 Yi 以前に事象が発生していない被験者 j の集合（被験者 i 自身を含む）に対するものである。明らかに、0 < Li(β) ≤ 1 である。これは部分尤度（偏尤度とも）（英語版）であり、時間経過によるハザードの変化をモデル化することなく、共変量の影響を推定することができる。 被験者を統計的に互いに独立であるかのように扱うと、実現したすべての事象の同時確率は、事象の発生を Ci = 1 で示すと、次のような部分尤度となる。 L ( β ) = ∏ i : C i = 1 L i ( β ) {\displaystyle L(\beta )=\prod _{i:C_{i}=1}L_{i}(\beta )} これに対応する対数部分尤度は、 ℓ ( β ) = ∑ i : C i = 1 ( X i ⋅ β − log ⁡ ∑ j : Y j ≥ Y i θ j ) {\displaystyle \ell (\beta )=\sum _{i:C_{i}=1}\left(X_{i}\cdot \beta -\log \sum _{j:Y_{j}\geq Y_{i}}\theta _{j}\right)} である。 この関数は β 上で最大化され、モデルパラメータの最大部分尤度推定量を生成することができる。 部分スコア関数（英語版）は、 ℓ ′ ( β ) = ∑ i : C i = 1 ( X i − ∑ j : Y j ≥ Y i θ j X j ∑ j : Y j ≥ Y i θ j ) {\displaystyle \ell ^{\prime }(\beta )=\sum _{i:C_{i}=1}\left(X_{i}-{\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}\right)} であり、部分対数尤度のヘッセ行列は、 ℓ ′ ′ ( β ) = − ∑ i : C i = 1 ( ∑ j : Y j ≥ Y i θ j X j X j ′ ∑ j : Y j ≥ Y i θ j − [ ∑ j : Y j ≥ Y i θ j X j ] [ ∑ j : Y j ≥ Y i θ j X j ′ ] [ ∑ j : Y j ≥ Y i θ j ] 2 ) {\displaystyle \ell ^{\prime \prime }(\beta )=-\sum _{i:C_{i}=1}\left({\frac {\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}X_{j}^{\prime }}{\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}}}-{\frac {\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}\right]\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}X_{j}^{\prime }\right]}{\left[\sum _{j:Y_{j}\geq Y_{i}}\theta _{j}\right]^{2}}}\right)} である。 このスコア関数とヘッセ行列を用いると、ニュートン＝ラフソンアルゴリズムを使用して部分尤度を最大化することができる。β の推定量で評価されるヘッセ行列の逆行列は、推定量の近似分散共分散行列として使用でき、回帰係数の近似標準誤差を生成するために使用できる。 時間データに同値がある場合に対処するために、いくつかの方法が提案されている。Breslow法（Breslow's method）は、同値が存在する場合でも、上述の手順を変更せずに使用するアプローチである。より良い結果が得られると考えられる別のアプローチとして、Efron法（Efron's method）がある。tj を一意の時間とし、 Hj を Yi = tj かつ Ci = 1となるインデックス i の集合とし、mj = |Hj| とする。Efronのアプローチは、次の部分尤度を最大化する。 L ( β ) = ∏ j ∏ i ∈ H j θ i ∏ ℓ = 0 m − 1 [ ∑ i : Y i ≥ t j θ i − ℓ m ∑ i ∈ H j θ i ] . {\displaystyle L(\beta )=\prod _{j}{\frac {\prod _{i\in H_{j}}\theta _{i}}{\prod _{\ell =0}^{m-1}\left[\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}\right]}}.} 対応する対数部分尤度は、 ℓ ( β ) = ∑ j ( ∑ i ∈ H j X i ⋅ β − ∑ ℓ = 0 m − 1 log ⁡ ( ∑ i : Y i ≥ t j θ i − ℓ m ∑ i ∈ H j θ i ) ) , {\displaystyle \ell (\beta )=\sum _{j}\left(\sum _{i\in H_{j}}X_{i}\cdot \beta -\sum _{\ell =0}^{m-1}\log \left(\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}\right)\right),} スコア関数は、 ℓ ′ ( β ) = ∑ j ( ∑ i ∈ H j X i − ∑ ℓ = 0 m − 1 ∑ i : Y i ≥ t j θ i X i − ℓ m ∑ i ∈ H j θ i X i ∑ i : Y i ≥ t j θ i − ℓ m ∑ i ∈ H j θ i ) , {\displaystyle \ell ^{\prime }(\beta )=\sum _{j}\left(\sum _{i\in H_{j}}X_{i}-\sum _{\ell =0}^{m-1}{\frac {\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}X_{i}}{\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}}}\right),} そしてヘッセ行列は、 ℓ ′ ′ ( β ) = − ∑ j ∑ ℓ = 0 m − 1 ( ∑ i : Y i ≥ t j θ i X i X i ′ − ℓ m ∑ i ∈ H j θ i X i X i ′ ϕ j , ℓ , m − Z j , ℓ , m Z j , ℓ , m ′ ϕ j , ℓ , m 2 ) , {\displaystyle \ell ^{\prime \prime }(\beta )=-\sum _{j}\sum _{\ell =0}^{m-1}\left({\frac {\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}X_{i}^{\prime }-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}X_{i}X_{i}^{\prime }}{\phi _{j,\ell ,m}}}-{\frac {Z_{j,\ell ,m}Z_{j,\ell ,m}^{\prime }}{\phi _{j,\ell ,m}^{2}}}\right),} であり、ここに ϕ j , ℓ , m = ∑ i : Y i ≥ t j θ i − ℓ m ∑ i ∈ H j θ i {\displaystyle \phi _{j,\ell ,m}=\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}} Z j , ℓ , m = ∑ i : Y i ≥ t j θ i X i − ℓ m ∑ i ∈ H j θ i X i . {\displaystyle Z_{j,\ell ,m}=\sum _{i:Y_{i}\geq t_{j}}\theta _{i}X_{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}X_{i}.} となる。なお、Hj が空の場合（時刻 tj のすべての観測値が打ち切られる）、これらの式の被加数はゼロとして扱われることに注意すること。 時間依存変数、時間依存層、および被験者ごとの複数の事象への拡張は、AndersenとGillの計数過程の定式化によって組み入れることができる。時変予測変数[訳語疑問点]を用いたハザードモデルの使用例として、失業保険の失業期間に対する効果の推定がある。 時変共変量（英語版）（すなわち予測変数）を許容することに加えて、Coxモデルは時変係数に一般化することができる。つまり、治療の比例効果は時間とともに変化する可能性がある。たとえば、ある薬剤は罹患後1ヶ月以内に投与すると非常に効果的だが、時間が経つにつれて効果が低下する場合がある。次に、係数が時間とともに変化しない（定常性）という仮説を検証することができる。詳細とソフトウェア（Rパッケージ）はMartinussen and Scheike (2006)から入手できる。信頼性計算では、時変共変量を使用したCoxモデルの応用が考えられている。 これに関連して、共変量の効果を加法ハザードで指定することが理論的に可能であることも言及されている。つまり、 λ ( t | X i ) = λ 0 ( t ) + β 1 X i 1 + ⋯ + β p X i p = λ 0 ( t ) + X i ⋅ β {\displaystyle \lambda (t|X_{i})=\lambda _{0}(t)+\beta _{1}X_{i1}+\cdots +\beta _{p}X_{ip}=\lambda _{0}(t)+X_{i}\cdot \beta } を指定することである。 （対数）尤度最大化を目的とした状況でこのような加法ハザードモデル（additive hazards models）を使用する場合は、 λ ( t ∣ X i ) {\displaystyle \lambda (t\mid X_{i})} を非負値に制限するように注意しなければならない。おそらくこの複雑さのために、このようなモデルはほとんど見られない。目的が最小二乗であれば、非負の制限は厳密には必要ない。 ベースライン・ハザードが特定の形に従うと仮定する理由がある場合、Coxモデルを特殊化することができる。この場合、ベースライン・ハザード λ 0 ( t ) {\displaystyle \lambda _{0}(t)} は与えられた関数で置き換えられる。たとえば、ハザード関数をワイブル・ハザード関数と仮定すると、「ワイブル比例ハザードモデル」（Weibull proportional hazards model）となる。 ちなみに、ワイブル・ベースライン・ハザードを使用することは、モデルが比例ハザードモデルと加速故障時間モデル（英語版）の両方を満足する唯一の状況である。 一般用語のパラメトリック比例ハザードモデル（parametric proportional hazards models）は、ハザード関数が指定されている比例ハザードモデルを説明するために使用できる。対照的に、Cox比例ハザードモデル（Cox proportional hazards model）はセミパラメトリック・モデル（英語版）と呼ばれることもある。 一部の著者は、基礎となるハザード関数を指定している場合でも、（パラメトリック比例ハザードモデルでなく）Cox比例ハザードモデルという用語を使用して、この分野全体がデイヴィッド・コックスに負っている恩義を認めている。 「Cox回帰モデル」（Cox regression model）（比例ハザードを除く）という用語は、時間依存因子を含むようにCoxモデルを拡張したものを表すために使われることがある。しかし、Cox比例ハザードモデルはそれ自体が回帰モデルとして記述できるため、この使い方は潜在的に曖昧である。 比例ハザードモデルとポアソン回帰（英語版）モデルの間には関係があり、ポアソン回帰のソフトウェアで近似比例ハザードモデルを適合させるために使用されることがある。これを行う一般的な理由は、計算がはるかに速くなるということである。このことは、コンピュータの速度が遅い時代にはより重要であったが、特に大きなデータセットや複雑な問題には依然として有益である。Laird and Olivier（1981）は、数学的な詳細を説明した。彼らは『我々は（ポアソンモデルが）真であると仮定しているのではなく、単に尤度を導出するための装置として使用するだけである』と述べている。McCullagh and Nelderの一般化線形モデルに関する本には、比例ハザードモデルから一般化線形モデルへの変換に関する章がある。 高次元では、共変量の数 p がサンプルサイズ n に比べて大きい場合、Lasso法は古典的なモデル選択戦略の1つである。Tibshirani（1997）は、比例ハザード回帰パラメータのLasso法を提案した。回帰パラメータ β のLasso推定量は、L1ノルム制約の下でのCox部分対数尤度の逆数の最小化として定義される。 ℓ ( β ) = ∑ j ( ∑ i ∈ H j X i ⋅ β − ∑ ℓ = 0 m − 1 log ⁡ ( ∑ i : Y i ≥ t j θ i − ℓ m ∑ i ∈ H j θ i ) ) + λ ‖ β ‖ 1 , {\displaystyle \ell (\beta )=\sum _{j}\left(\sum _{i\in H_{j}}X_{i}\cdot \beta -\sum _{\ell =0}^{m-1}\log \left(\sum _{i:Y_{i}\geq t_{j}}\theta _{i}-{\frac {\ell }{m}}\sum _{i\in H_{j}}\theta _{i}\right)\right)+\lambda \|\beta \|_{1},} 最近、このトピックについて理論的な進歩があった。 ポータル 数学 加速故障時間モデル（英語版） 10％ルール（英語版） ワイブル分布 尤度方程式