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)Name,Value ペア引数で指定された追加オプションで、係数推定値のベクトルを返します。
例
標本データを読み込みます。
load('lightbulb.mat');電球データの 1 列目には 2 種類の電球の寿命 (時間単位) が含まれています。2 列目には、電球が蛍光灯か白熱灯であるかを示すバイナリ変数が含まれています。0 は電球が蛍光灯であることを、1 は白熱灯であることを示します。3 列目には打ち切り情報が含まれます。ここで、0 は電球が故障するまで観測されたことを示し、1 は電球が打ち切られたことを示します。
電球の寿命のコックス比例ハザード モデルを当てはめ、打ち切りも考慮します。予測子変数は電球のタイプです。
b = coxphfit(lightbulb(:,2),lightbulb(:,1), ... 'Censoring',lightbulb(:,3))
b = 4.7262
ハザード率の推定値は = 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: [100×1 double]
devres: [100×1 double]
martres: [100×1 double]
schres: [100×1 double]
sschres: [100×1 double]
scores: [100×1 double]
sscores: [100×1 double]
LikelihoodRatioTestP: 6.6613e-16
係数推定値の共分散行列 covb に含まれる 1 つの値は、この例の係数推定値の分散と等しくなります。係数推定値 beta は b と同じで、0.9409 になります。係数推定値の標準誤差 se は 0.1256 で、分散 0.0158 の平方根になります。 統計量 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')
近似モデルによって、実際の分布の生存時間関数に近い推定値が示されます。
入力引数
名前と値の引数
出力引数
詳細
アルゴリズム
階層のベースライン累積ハザード率 (
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.
