Main Content

coxphfit

コックス比例ハザード回帰

ページ内をすべて折りたたむ

構文

b = coxphfit(X,T)
b = coxphfit(X,T,Name,Value)
[b,logl,H,stats] = coxphfit(___)

説明

例

b = coxphfit(X,T) は、予測子 X に対する、観測された応答 T のコックス比例ハザード回帰の係数推定値が格納されている、p 行 1 列のベクトル b を返します。T は n 行 1 列のベクトルまたは n 行 2 列の行列、X は n 行 p 列の行列です。

このモデルには定数項がありません。また、X には 1 で構成される列を含めることはできません。

例

b = coxphfit(X,T,Name,Value) は、1 つ以上の Name,Value ペア引数で指定された追加オプションで、係数推定値のベクトルを返します。

例

また、[b,logl,H,stats] = coxphfit(___) は、対数尤度 logl、追加の統計を含む stats 構造体、最初の列に T 値、2 番目の列に推定ベースライン累積ハザードを含む 2 列の行列 H を返します。前の構文の入力引数のいずれかを使用できます。

例

すべて折りたたむ

ライブ スクリプトを開く

標本データを読み込みます。 

load('lightbulb.mat');

電球データの 1 列目には 2 種類の電球の寿命 (時間単位) が含まれています。2 列目には、電球が蛍光灯か白熱灯であるかを示すバイナリ変数が含まれています。0 は電球が蛍光灯であることを、1 は白熱灯であることを示します。3 列目には打ち切り情報が含まれます。ここで、0 は電球が故障するまで観測されたことを示し、1 は電球が打ち切られたことを示します。

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

b = coxphfit(lightbulb(:,2),lightbulb(:,1), ...
'Censoring',lightbulb(:,3))
b = 4.7262

ハザード率の推定値は eb = 112.8646 です。つまり、白熱灯電球のハザードは、蛍光灯電球のハザードの 112.86 倍であることを意味します。

ライブ スクリプトを開く

標本データを読み込みます。 

load('lightbulb.mat');

データの 1 列目には 2 種類の電球の寿命 (時間単位) が含まれています。2 列目には、電球が蛍光灯か白熱灯であるかを示すバイナリ変数が含まれています。1 は電球が蛍光灯であることを、0 は白熱灯であることを示します。3 列目には打ち切り情報が含まれます。ここで、0 は電球が故障するまで観測されることを示し、1 はアイテム (電球) が打ち切られることを示します。

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

b = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3))
b = 4.7262

coxphfit で係数の推定に使用するアルゴリズムの既定の制御パラメーターを表示します。 

statset('coxphfit')
ans = struct with fields:
          Display: 'off'
      MaxFunEvals: 200
          MaxIter: 100
           TolBnd: []
           TolFun: 1.0000e-08
       TolTypeFun: []
             TolX: 1.0000e-08
         TolTypeX: []
          GradObj: []
         Jacobian: []
        DerivStep: []
      FunValCheck: []
           Robust: []
     RobustWgtFun: []
           WgtFun: []
             Tune: []
      UseParallel: []
    UseSubstreams: []
          Streams: {}
        OutputFcn: []

オプションを異なる名前で保存し、結果の表示方法と最大反復回数 (Display および MaxIter) を変更します。 

coxphopt = statset('coxphfit');
coxphopt.Display = 'final';
coxphopt.MaxIter = 50;

新しいアルゴリズム パラメーターを使用して coxphfit を実行します。 

b = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3),'Options',coxphopt)
Successful convergence: Norm of gradient less than OPTIONS.TolFun

b = 4.7262

coxphfit は、最後の反復に関するレポートを表示します。最大反復回数を変更しても、係数推定値には影響しません。

ライブ スクリプトを開く

予測子 X に依存するワイブル データを生成します。 

rng('default') % for reproducibility
X = 4*rand(100,1);
A = 50*exp(-0.5*X); 
B = 2;
y = wblrnd(A,B);

応答値は、予測子変数 X と形状パラメーター 2 によって変化するスケール パラメーターをもつワイブル分布から生成されます。

コックス比例ハザード モデルを当てはめます。 

[b,logL,H,stats] = coxphfit(X,y);
[b logL]
ans = 1×2

    0.9409 -331.1479

係数推定値は 0.9409 で、対数尤度値 は –331.1479 です。

モデルの統計量を要求します。 

stats
stats = struct with fields:
                    covb: 0.0158
                    beta: 0.9409
                      se: 0.1256
                       z: 7.4889
                       p: 6.9462e-14
                   csres: [100x1 double]
                  devres: [100x1 double]
                 martres: [100x1 double]
                  schres: [100x1 double]
                 sschres: [100x1 double]
                  scores: [100x1 double]
                 sscores: [100x1 double]
    LikelihoodRatioTestP: 6.6613e-16

係数推定値の共分散行列 covb に含まれる 1 つの値は、この例の係数推定値の分散と等しくなります。係数推定値 beta は b と同じで、0.9409 になります。係数推定値の標準誤差 se は 0.1256 で、分散 0.0158 の平方根になります。z 統計量 z は、beta/se = 0.9409/0.1256 = 7.4880 です。p 値 p は、X の効果が有意であることを示します。

ベースライン生存時間関数のコックス推定値を既知のワイブル関数と共にプロットします。 

stairs(H(:,1),exp(-H(:,2)),'LineWidth',2)
xx = linspace(0,100);
line(xx,1-wblcdf(xx,50*exp(-0.5*mean(X)),B),'color','r','LineWidth',2)
xlim([0,50])
legend('Estimated Survivor Function','Weibull Survivor Function')

Figure contains an axes object. The axes object contains 2 objects of type stair, line. These objects represent Estimated Survivor Function, Weibull Survivor Function.

近似モデルによって、実際の分布の生存時間関数に近い推定値が示されます。

入力引数

すべて折りたたむ

予測子変数に関する観測値。n 個の観測値のそれぞれについて p 個の予測子の n 行 p 列の行列として指定します。

このモデルには定数項がないため、X には 1 で構成される列を含めることはできません。

X、T、あるいは 'Frequency' または 'Strata' の値に NaN 値が含まれている場合、coxphfit はコックス モデルを当てはめるときに、NaN 値がある行をすべてのデータから削除します。

データ型: double

イベントが発生するまでの時間のデータ。n 行 1 列のベクトルまたは 2 列の行列を指定します。

  • T が n 行 1 列のベクトルの場合、イベントが発生するまでの時間のデータを右側打ち切りしたイベント時間を表します。

  • T が n 行 2 列の行列の場合、各行は時間依存共変量のリスク区間 (start,stop] を計数過程形式で表します。1 列目は開始時間、2 列目は停止時間です。たとえば、共変量が時間に依存するコックス比例ハザード モデルを参照してください。

X、T、あるいは 'Frequency' または 'Strata' の値に NaN 値が含まれている場合、coxphfit はコックス モデルを当てはめるときに、NaN 値がある行をすべてのデータから削除します。

データ型: single | double

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: 'Baseline',0,'Censoring',censoreddata,'Frequency',freq の指定により、coxphfit では、ベクトル censoreddata の打ち切り情報と、ベクトル freq で与えられる T および X での観測の頻度を考慮して、0 を基準にしたベースライン ハザード率を計算します。

係数の初期値。'B0' と数値ベクトルから構成されるコンマ区切りの値として指定します。

データ型: double

ベースライン ハザードの計算対象となる X 値。'Baseline' とスカラー値で構成されるコンマ区切りのペアとして指定します。

既定値は mean(X) であるため、X のハザード率は h(t)*exp((X-mean(X))*b) です。0 を基準にしてベースラインを計算するには「0」を入力します。これにより、X のハザード率は h(t)*exp(X*b) になります。ベースラインを変更しても係数推定値には影響しませんが、ハザード率が変わります。

例: 'Baseline',0

データ型: double

打ち切りのインジケーター。'Censoring' と、T と同じサイズの論理値配列で構成されるコンマ区切りのペアで指定します。右側打ち切りの観測値の場合は 1、完全に観測された観測値の場合は 0 を使用します。既定の設定では、すべての観測値が完全に観測されます。たとえば、打ち切りデータのコックス比例ハザード モデルを参照してください。

例: 'Censoring',cens

データ型: logical

観測値の頻度または重み。'Frequency' と配列から構成されるコンマ区切りのペアとして指定します。この配列のサイズは、非負のスカラー値が格納されている T のサイズと同じです。この配列には、観測値の頻度に対応する整数値、または観測値の重みに対応する非負値を格納できます。

X、T、あるいは 'Frequency' または 'Strata' の値に NaN 値が含まれている場合、coxphfit はコックス モデルを当てはめるときに、NaN 値がある行をすべてのデータから削除します。

既定は、X および T の各行に対して 1 です。

例: 'Frequency',w

データ型: double

階層化変数。実数値の行列から構成されるコンマ区切りのペアとして指定します。この行列の行数は T と同じでなければならず、各行は 1 つの観測値に対応しなければなりません。

X、T、あるいは 'Frequency' または 'Strata' の値に NaN 値が含まれている場合、coxphfit はコックス モデルを当てはめるときに、NaN 値がある行をすべてのデータから削除します。

既定の [] は階層化変数なしです。

例: 'Strata',Gender

データ型: single | double

同順位の故障時間の処理方法。'Ties' と 'breslow' (Breslow の方法) または 'efron' (Efron の方法) のいずれかから構成されるコンマ区切りのペアとして指定します。

例: 'Ties','efron'

b を推定するために使用される反復アルゴリズムのアルゴリズム制御パラメーター。'Options' と構造体で構成されるコンマ区切りのペアとして指定します。statset を呼び出すと、この引数が作成されます。パラメーター名と既定値には、「statset('coxphfit')」と入力します。オプションを新しい名前で設定し、その名前を名前と値のペア引数で使用できます。

例: 'Options',statset('coxphfit')

出力引数

すべて折りたたむ

コックス比例ハザード回帰の係数推定値。p 行 1 列のベクトルとして返されます。

近似されたモデルの対数尤度。スカラーとして返されます。

対数尤度値を使用してさまざまなモデルを比較し、モデルの項の効果の有意性を評価できます。

T の値で評価された推定ベースライン累積ハザード率。以下のいずれかとして返されます。

  • モデルが階層化されていない場合、H は 2 列の行列です。行列の最初の列には T 値が含まれ、2 列目には累積ハザード率の推定値が含まれます。

  • モデルが階層化されている場合、H は (2+k) 列の行列です。最後の k 列は、名前と値のペアの引数 Strata を使用した階層化変数に対応します。

係数統計量。以下のフィールドを含む構造体として返されます。

beta係数推定値 (b と同じ)
se係数推定値 b の標準誤差
zb の z 統計量 (b を標準誤差で除算したもの)
pb の p 値
covb

b に対して推定された共分散行列

csres

Cox-Snell 残差

devres逸脱度残差
martresマルチンゲール残差
schresシェーンフィールド残差
sschresスケーリングされたシェーンフィールド残差
scoresスコア残差
sscoresスケーリングされたスコア残差

coxphfit は、Cox-Snell 残差、マルチンゲール残差および逸脱度残差を、各行に観測値がある列ベクトルとして返します。シェーンフィールド残差、スケーリングされたシェーンフィールド残差、スコア残差およびスケーリングされたスコア残差は、X と同じサイズの行列として返されます。打ち切られたデータのシェーンフィールド残差およびスケーリングされたシェーンフィールド残差は NaN になります。

詳細

すべて折りたたむ

コックス比例ハザード回帰

コックス比例ハザード回帰は、生存率の推定値を調整して交絡変数の影響を除去したり、予測子変数の影響を定量化するためのセミパラメトリック手法です。この手法は、注釈変数および交絡変数の影響を一般的なベースライン ハザード関数 h0(t) の乗数として表します。

0 に対するベースラインの場合、このモデルは次に対応します。

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

ここで、Xi=(xi1,xi2,⋯,xip) は i 番目の被験者についての予測子変数、h(Xi,t) は時間 t における Xi についてのハザード率、h0(t) はベースライン ハザード率関数です。ベースライン ハザード関数はコックス比例ハザード回帰関数のノンパラメトリック部分ですが、予測子変数の影響は対数線形回帰です。このモデルは、ベースライン ハザード関数が時間 t に依存しているが、予測子変数は時間に依存していないという仮定に基づいています。階層化および時間依存変数に向けた拡張、同時発生イベント、観測値の重みなどの詳細については、コックス比例ハザード モデルを参照してください。

アルゴリズム

  • 階層のベースライン累積ハザード率 (H) を計算する場合は、階層の入力データに完全に観測された観測値が少なくとも 1 つは含まれていなければなりません。階層に打ち切られた観測値しか含まれていない場合、出力 H の最初の 2 列に NaN が格納され、残りの列に階層情報が格納されます。

    R2022a より前: 階層に打ち切られた観測値しか含まれていない場合、H にゼロの行が格納され、階層情報は格納されません。

参照

[1] Cox, D.R., and D. Oakes. Analysis of Survival Data. London: Chapman & Hall, 1984.

[2] Lawless, J. F. Statistical Models and Methods for Lifetime Data. Hoboken, NJ: Wiley-Interscience, 2002.

[3] Kleinbaum, D. G., and M. Klein. Survival Analysis. Statistics for Biology and Health. 2nd edition. Springer, 2005.

拡張機能

バージョン履歴

R2006a より前に導入

すべて展開する

参考

| |

トピック