Main Content

CoxModel

コックス比例ハザード モデル

R2021a 以降

説明

コックス比例ハザード モデルは、寿命や故障時間のデータに関連するモデルです。基本的なコックス モデルには、ハザード関数 h0(t) とモデル係数 b が含まれており、予測子 X についての時間 t におけるハザード率は次のようになります。

h(Xi,t)=h0(t)exp[j=1pxijbj],

係数 b は時間に依存しません。作成関数 fitcox は、モデル係数 b とハザード率 h0(t) の両方を推定し、それらを結果の CoxModel オブジェクトにプロパティとして格納します。

完全なコックス モデルは、基本的なモデルを拡張したもので、異なるベースラインを基準とするハザードや階層化変数などが含まれています。コックス比例ハザード モデルの拡張を参照してください。

作成

CoxModel オブジェクトを作成するには、fitcox を使用します。

プロパティ

すべて展開する

モデルを当てはめるときに指定されたベースライン ハザード。実数スカラーとして指定されます。コックス モデルは相対ハザード モデルであるため、与えられたデータのハザードを比較する基準となるベースラインが必要です。既定値は mean(X) (階層化モデルの場合は各階層化における平均) であるため、X のハザード率は h(t)*exp((X-mean(X))*b) です。0 を基準にしてベースラインを計算するには「0」を入力します。これにより、X のハザード率は h(t)*exp(X*b) になります。ベースラインを変更しても係数には影響しません。

データ型: double

係数推定値の共分散行列。行数が予測子の数と等しい正方行列として指定されます。

データ型: double

係数と関連する統計。次の 4 列の table として指定されます。

  • Beta — 係数推定値

  • SE — 係数推定値の標準誤差

  • zStat — z 統計量

  • pValue — 係数の p 値 (ゼロ Beta と比較)

table の各行は 1 つの予測子に対応します。

ベクトルとしてこれらの列のいずれかを取得するには、ドット表記を使ってプロパティにインデックスを付けます。たとえば、coxMdl オブジェクトの推定される係数ベクトルは次のとおりです。

beta = coxMdl.Coefficients.Beta

係数に対する他の検定を実行するには、linhyptest を使用します。

データ型: table

当てはめに使用したモデルの表現。ウィルキンソンの表記法による式として指定されます。ウィルキンソンの表記法を参照してください。たとえば、複数の予測子を含めるには次を使用します。

'X ~ a + b + … + c'

ここで、変数 abc のそれぞれが 1 つの予測子を表します。これらの変数は table X の列名です。

推定ベースライン累積ハザード。double の行列として指定されます。累積ハザードは学習で定義された時間点において評価されます。

Hazard には少なくとも 2 つの列があります。1 列目には時間値が格納され、残りの列にはリストされた各時間における累積ハザードが格納されます。

  • 非階層化モデルの場合、Hazard は 2 列です。

  • 階層化モデルの場合、Hazard に階層化レベルの一意の組み合わせごとに追加の列が含まれます。それぞれの階層化に対応する Hazard(:,1) の一意の時間値が Hazard(:,2) ではエントリ 0 で区切られます。階層化モデルは、名前と値の引数 'Stratification' を使用して学習させたモデルです。

理論的には、時間 t における累積ハザードは –log(1 – cdf(t)) です。経験的累積ハザードは次のようになります。

H0^(t)=tith0^(ti)=tit1jRiexp(β·xj),

ここで、Ri は時間 ti におけるリスク集合、つまり故障のリスクがある観測値です。部分尤度関数を参照してください。

データ型: double

モデルがヌル モデルに対して有意であるかどうかを示す P 値。実数スカラーとして指定されます。

このプロパティには、ヌル モデル、つまりすべての係数が 0 に等しいモデルに対して実行した尤度比検定の p 値が格納されます。

この尤度比検定では、係数推定値におけるデータとすべての係数が 0 である場合のデータで尤度関数を比較します。この比較により、学習させたモデルがすべての係数が 0 に等しいモデルに対して有意であるかどうかを判別するのに使用できる検定統計量が得られます。帰無仮説はヌル モデルと学習させたモデルに差異がないことであるため、p 値が有意であれば、学習させたモデルが有意であることを意味します。

データ型: double

係数推定値における尤度関数の対数。実数スカラーとして指定されます。

データ型: double

モデルの予測子 (係数) の数。正の整数として指定されます。

データ型: double

モデルの当てはめに使用された予測子の名前。文字ベクトルの cell 配列として指定されます。モデルを table のデータで学習させた場合、予測子の名前は table の列の名前になります。それ以外の場合は、予測子の名前は X1X2 のようになります。

データ型: cell

各共変量が比例ハザードの仮定を満たすかどうかを示す P 値。予測子ごとに 1 つのエントリをもつ実数ベクトルとして指定されます。

コックス モデルは比例ハザードの仮定に依存し、つまり、任意の 2 つのデータ点 X1 および X2 について、hazard(X1)/hazard(X2) が一定になります。予測子が時間に依存する場合、この仮定に違反することがあります。たとえば、予測子が年齢に対応する場合、一般に年齢が増すほどハザードが高まります。

この仮定の検定では、スケーリングされたシェーンフィールド残差を使用します。これは、Grambsch と Therneau が[1]で導き出したものです。

帰無仮説は、各係数が比例ハザードの仮定を満たすことです。p 値が有意であれば、特定の係数が比例ハザードの仮定に違反していることを意味します。検定は各係数について実行されるため、このプロパティは係数と同じ数の要素をもつベクトルになります。

データ型: double

モデル全体が比例ハザードの仮定を満たすかどうかを示す P 値。実数スカラーとして指定されます。

コックス モデルは比例ハザードの仮定に依存し、つまり、任意の 2 つのデータ点 X1 および X2 について、hazard(X1)/hazard(X2) が一定になります。予測子が時間に依存する場合、この仮定に違反することがあります。たとえば、予測子が年齢に対応する場合、一般に年齢が増すほどハザードが高まります。

この仮定の検定では、スケーリングされたシェーンフィールド残差を使用します。これは、Grambsch と Therneau が[1]で導き出したものです。

帰無仮説は、モデル全体が比例ハザードの仮定を満たすことです。p 値が有意であれば、モデル全体が比例ハザードの仮定を満たしていないことを意味します。

データ型: double

各種の残差。各列がそれぞれの残差に対応する 7 列の table として指定されます。

  • 'CoxSnell' — 観測値 X(i) の Cox-Snell 残差は、時間 i における累積ハザード (cumHazard(i)) に X(i) のハザードを乗算して csres(i) = cumHazard(i) * exp(X(i) * Beta) と定義されます。Beta は、Coefficients に格納されている当てはめた Beta ベクトルです。

  • 'Deviance' — 逸脱度残差は、マルチンゲール残差を使用して D(i) = sign(M(i))*sqrt(–2*[M(i) + delta(i)*log(delta(i)–M(i)))] と定義されます。ここで、D(i)i 番目の逸脱度残差、M(i)i 番目のマルチンゲール残差であり、delta(i) はデータ点 i が死亡したかどうかを示します。

  • 'Martingale' — 点 X(i) のマルチンゲール残差は、delta(i) – CoxSnell(i) となります。ここで、delta(i)X(i) が死亡したかどうかを示し、CoxSnell(i)i における Cox-Snell 残差です。マルチンゲール残差は、X(i) の真の死亡数からモデルに基づく予測死亡数を減算した差と見なすことができます。

  • 'Schoenfeld' — シェーンフィールド残差は、scres(i,j) = X(i,j) – M(Beta,i,j) と定義されます。ここで、X(i,j) は観測値 ij 番目の要素であり、M(Beta,i,j) は、時間 i に残された生存観測値数を仮定した、X(i,j) の期待値です。シェーンフィールド残差は、時間 i での真の死亡観測値と、残りの生存観測値を仮定した、時間 i での確認にモデルが期待する死亡観測値の差と見なすことができます。各共変量について残差が計算されるため、学習したパラメーターと同じ数の列が含まれます。この残差は、死亡が観測された時間と観測値に対してのみ有効です。打ち切られた観測値についてはすべて、対応する残差が NaN になります。

  • 'ScaledSchoenfeld' — スケーリングされたシェーンフィールド残差は、学習した係数の分散でシェーンフィールド残差をスケーリングしたものです。シェーンフィールド残差と同様に、スケーリングされた残差は死亡に至っていない観測値と時間については定義されず、そのような点の残差は NaN になります。

  • 'Score' — スコア残差は、scores(i,t) = integral_{0}^{t}(X(i,u) – Xbar(u)) dScres(i,u) と定義されます。ここで、Schres(i) は観測値 i のシェーンフィールド残差であり、Xbar は時間 u に生存している観測値の平均です。各共変量について残差が計算されるため、学習したパラメーターと同じ数の列が含まれます。

  • 'ScaledScore' — スケーリングされたスコア残差は、当てはめた係数の共分散でスコア残差をスケーリングしたものです。

Residuals には学習データと同じ数の行が含まれます。

データ型: table

応答変数名。文字ベクトルを指定します。応答値が table に格納されたモデルの場合、応答変数名は関連する table の列の名前になります。それ以外の場合は、ResponseName'T' です。

データ型: char

係数推定値の標準誤差。実数ベクトルとして指定されます。StandardError は、行列 CoefficientCovariance の対角の平方根になります。

データ型: double

学習における入力階層化の一意の組み合わせの配列。次のいずれかのデータ型として指定されます。

  • 数値配列 — すべての階層化変数が数値。

  • string 配列 — すべての変数が string。

  • string の cell 配列 — すべての変数が cell string。

  • categorical 配列 — すべての変数が categorical。

  • cell 配列 — 変数の型が混在。

XT のいくつかのデータについて、次の表に Stratification の内容の例を示します。

入力データ結果として得られる階層化
double mdl = fitcox(X,T,'Stratification',[1 2; 1 2; 2 2; 2 2]); [1 2; 2 2]
stringmdl = fitcox(X,T,'Stratification',["cat";"dog";"dog";"bird"]);["cat"; "dog"; "bird"]
cell stringmdl = fitcox(X,T,'Stratification',{'cat';'dog';'dog';'bird'});{'cat';'dog';'bird'}
categoricalmdl = fitcox(X,T,'Stratification',categorical([1;2;3;4]));categorical([1;2;3;4])
混在する型

data = table(X,T,[1;2;1;3],{'cat';'cat';'dog';'dog'},'VariableNames',{'X','T','S1','S2');

mdl = fitcox(data,'T','Stratification',{'S1','S2'});

{1 'cat'; 2 'cat'; 1 'dog'; 3 'dog'}

データ型: double | char | string | cell | categorical

データの当てはめに関する情報。次の 4 列の table として指定されます。

  • Class — 予測子のクラス。

  • Range — 予測子がカテゴリカルでない場合は最小値と最大値、予測子がカテゴリカルの場合はすべてのカテゴリのリスト。

  • InModel — 予測子がモデルで使用されているかどうかを示す logical。応答変数はモデルに含まれません。学習に使用された予測子変数はモデルに含まれます。

  • IsCategorical — 予測子が学習においてカテゴリカルとして扱われたかどうかを示す logical。

モデルにカテゴリカル予測子がなく、モデルの当てはめに式が使用されなかった場合は、VariableInfo の行数はモデルの予測子の数になります。それ以外の場合は、行数は PredictorNames 内の要素の数と同じになります。

データ型: table

オブジェクト関数

coefciコックス比例ハザード モデルの係数の信頼区間
discardResidualsコックス モデルから残差を削除
hazardratioコックス モデルのベースラインに対するハザードの推定
linhyptestコックス モデルの係数に対する線形仮説検定
plotSurvivalコックス比例ハザード モデルの生存時間関数のプロット
survivalコックス比例ハザード モデルの生存の計算

すべて折りたたむ

ワイブル確率変数では、形状パラメーターが同じであればハザード率は比例します。ワイブル分布を参照してください。スケール パラメーターを a、形状パラメーターを b とすると、時間 t におけるハザード率は次のようになります。

babtb-1.

ワイブル分布から疑似無作為標本を生成します。スケール パラメーターは 1、5、および 1/3 で、形状パラメーターはいずれも B とします。

rng default % For reproducibility
B = 2;
A = ones(100,1);
data1 = wblrnd(A,B);
A2 = 5*A;
data2 = wblrnd(A2,B);
A3 = A/3;
data3 = wblrnd(A3,B);

データの table を作成します。予測子は 3 つの変数の型で、1、2、または 3 です。

predictors = categorical([A;2*A;3*A]);
data = table(predictors,[data1;data2;data3],'VariableNames',["Predictors" "Times"]);

データにコックス回帰を当てはめます。

mdl = fitcox(data,"Times")
mdl = 
Cox Proportional Hazards regression model

                     Beta        SE        zStat       pValue  
                    _______    _______    _______    __________

    Predictors_2    -3.5834    0.33187    -10.798    3.5299e-27
    Predictors_3     2.1668    0.20802     10.416    2.0899e-25


Log-likelihood: -1197.917

rates = exp(mdl.Coefficients.Beta)
rates = 2×1

    0.0278
    8.7301

電球の寿命のシミュレーション データが格納された lightbulb データ セットでコックス比例ハザード回帰を実行します。電球データの 1 列目には 2 種類の電球の寿命 (時間単位) が含まれています。2 列目には電球が蛍光灯と白熱灯のどちらであるかを示すバイナリ変数が含まれています。0 は電球が蛍光灯であることを示し、1 は白熱灯であることを示します。3 列目には打ち切り情報が含まれています。0 は電球が故障するまで観測されたことを示し、1 は観測が打ち切られたことを示します。

lightbulb データ セットを読み込みます。

load lightbulb

電球の寿命について、打ち切りを考慮してコックス比例ハザード モデルを当てはめます。予測子変数は電球のタイプです。

coxMdl = fitcox(lightbulb(:,2),lightbulb(:,1), ...
    'Censoring',lightbulb(:,3))
coxMdl = 
Cox Proportional Hazards regression model

           Beta       SE      zStat       pValue  
          ______    ______    ______    __________

    X1    4.7262    1.0372    4.5568    5.1936e-06


Log-likelihood: -212.638

exp(Beta) を評価して、蛍光灯電球と比較した白熱灯電球のハザード率を求めます。

hr = exp(coxMdl.Coefficients.Beta)
hr = 112.8646

ハザード率の推定値は eBeta = 112.8646 です。つまり、白熱灯電球の推定ハザードは蛍光灯電球のハザードの 112.86 倍であることを意味します。coxMdl.Coefficients.pValue の値が小さいことから、2 種類の電球のハザード率が同じになる、つまり Beta = 0 になる可能性は無視できます。

参照

[1] Grambsch, Patricia M., and Terry M. Therneau. Proportional Hazards Tests and Diagnostics Based on Weighted Residuals. Biometrika, vol. 81, no. 3, 1994, pp. 515–526. JSTOR, https://www.jstor.org/stable/2337123.

バージョン履歴

R2021a で導入