共変量が時間に依存するコックス比例ハザードモデル

この例では、生存データを計数過程形式に変換してから、共変量が時間に依存するコックス比例ハザードモデルを作成する方法を示します。

手順 1. 標準レイアウトと計数過程形式を比較する。

共変量が時間に依存するコックスモデルでは、生存データが標準レイアウトではなく計数過程形式になっている必要があります。標準レイアウトの生存データと計数過程形式の生存データの違いを調べるため、次の標本データを読み込みます。

load simplesurvivaldata

この標本データには、relapseS および relapseCP という 2 つのテーブルが含まれています。この 2 つのテーブルはそれぞれ、同じシンプルな生存データを標準レイアウトと計数過程形式で表しています。

標準レイアウトのデータを表示します。

relapseS

relapseS=2×5 table
    ID    Time    Censoring    Age    StopTreatment
    __    ____    _________    ___    _____________

    1       5         0        20          NaN     
    2      20         1        30           12

このデータは、時間の経過に従って治療状況が変化する 2 人の患者を表しています。患者 1 は、0 ～ 5 週目の期間に治療を受けておらず、この期間の終わりに再発しました。患者 2 は、0 ～ 12 週目の期間に治療を受けましたが、12 ～ 20 週目の期間には治療を受けませんでした。患者 2 は再発しておらず、20 週目の後で調査を打ち切りました。

次に、同じデータを計数過程形式で表示します。

relapseCP

relapseCP=3×6 table
    ID    tStart    tStop    Censoring    Age    TreatmentStatus
    __    ______    _____    _________    ___    _______________

    1        0        5          0        20            0       
    2        0       12          1        30            1       
    2       12       20          1        30            0

計数過程形式では、各イベント時間について各行が単一の値ではなくリスク区間 (tStart,tStop] を表しています。Censoring は、リスク区間の終わりでイベントが観測された場合は 0、観測されなかった場合は 1 です。TreatmentStatus は、時間に依存する共変量に対応しており、標準レイアウトの StopTreatment と同じ情報を表します。コックスモデルでは時間に依存する共変量が各リスク区間内で一定であると仮定することに注意してください。

手順 2. 標本データを読み込む。

次に、変換する標本データを読み込みます。

load survivaldatacp

この標本データにはテーブル labS が含まれ、各患者の反復測定を含むシミュレートされた生存データが標準レイアウトで格納されています。

標準レイアウトのシミュレートされた生存データを表示します。

labS

labS=6×7 table
    ID    Time    Censoring    Sex    Lab_0    Lab_50    Lab_100
    __    ____    _________    ___    _____    ______    _______

    1      46         0         1      0.3       NaN       NaN  
    2     138         1         0      0.2      0.23      0.39  
    3      94         0         1     0.18      0.22       NaN  
    4      50         0         0     0.21       0.2       NaN  
    5     106         0         0     0.25      0.21      0.42  
    6      98         0         0     0.21      0.22       NaN

標準レイアウトでは、テーブルの各行は 1 人の患者に関する情報を示します。

ID は患者の ID を示します。コックスモデルの入力に ID は含めません。ID は、データセットが計数過程形式に正しく変換されたことを確認するためにデータセットに含めます。
Time は、イベントが発生するまでの日数を表しており、応答変数に対応します。
Censoring には各患者の打ち切り情報が含まれています。1 は打ち切られたデータを示し、0 はイベントが発生するまでの正確な時間が観測期間の終わりで観測されたことを示します。
Sex は時間に依存しない予測子であり、1 は女性、0 は男性を示します。
Lab_0、Lab_50 および Lab_100 は、0、50 および 100 日目に測定された 3 つの連続する実験結果であり、時間に依存する予測子に対応します。

手順 3. 生存データを計数過程形式に変換する。

生存データ labS を計数過程形式に変換するため、以下のコードを実行します。このコードは、Time をリスク区間 (tStart,tStop] に変換し、時間に依存する予測子の 3 つのベクトル Lab_0、Lab_50 および Lab_100 を 1 つのベクトル Lab に結合します。

mTime = [0 50 100]; % Measurement time
threeLabs = [labS.Lab_0 labS.Lab_50 labS.Lab_100]; 
nLabMeasure = sum(sum(~isnan(threeLabs))); % Number of lab measurements
data = zeros(nLabMeasure,6); % One row for each observation
oID = 0; % Observation ID 
for i = 1 : size(labS,1)
    idx = find(mTime <= labS.Time(i));
    for j = 1 : length(idx)-1
        oID = oID + 1;
        data(oID,:) = [labS.ID(i) mTime(j:j+1) 1 labS.Sex(i) threeLabs(i,j)];
    end
    oID = oID + 1;
    data(oID,:) = [labS.ID(i) mTime(length(idx)) labS.Time(i) ... 
            labS.Censoring(i) labS.Sex(i) threeLabs(i,length(idx))];
end
labCP = table(data(:,1),data(:,2),data(:,3),data(:,4),data(:,5),data(:,6), ...
    'VariableNames', {'ID','tStart','tStop','Censoring','Sex','Lab'});

計数過程形式の生存データを表示します。

labCP

labCP=13×6 table
    ID    tStart    tStop    Censoring    Sex    Lab 
    __    ______    _____    _________    ___    ____

    1        0        46         0         1      0.3
    2        0        50         1         0      0.2
    2       50       100         1         0     0.23
    2      100       138         1         0     0.39
    3        0        50         1         1     0.18
    3       50        94         0         1     0.22
    4        0        50         1         0     0.21
    4       50        50         0         0      0.2
    5        0        50         1         0     0.25
    5       50       100         1         0     0.21
    5      100       106         0         0     0.42
    6        0        50         1         0     0.21
    6       50        98         0         0     0.22

計数過程形式では、テーブル labCP の各行は、1 つのリスク区間に対応する 1 つの観測値の情報を示します。コックスモデルではリスク区間 (tStart,tStop] で Lab が一定であると仮定することに注意してください。Censoring の値は、リスク区間の終わりでイベントが観測された場合は 0、観測されなかった場合は 1 です。

たとえば、患者 3 は 0 日目と 50 日目に 2 つの実験測定値があるので、計数過程形式では患者 3 に 2 行のデータがあります。コックスモデルでは、実験結果 0.18 および 0.22 がそれぞれ区間 (0,50] および (50,94] で一定であると仮定します。患者 3 の正確なイベント時間は 94 日目に観測されたので、Censoring は (0,50] では 1、(50,94] では 0 です。

手順 4. 長さがゼロのリスク区間を調整する。

リスク区間の長さがゼロである患者を探します。

idxInvalid = labCP.ID(find(labCP.tStart == labCP.tStop))

idxInvalid = 4

患者 4 のデータを確認します。

labCP(find(labCP.ID==idxInvalid),:)

ans=2×6 table
    ID    tStart    tStop    Censoring    Sex    Lab 
    __    ______    _____    _________    ___    ____

    4        0       50          1         0     0.21
    4       50       50          0         0      0.2

患者 4 のイベントは、50 日目の測定値で同時に発生しています。しかし、コックスモデルでは長さがゼロの区間が認められないので、(50,50] は無効なリスク区間です。リスク区間が有効になるように調整します。時間単位より小さい任意の値を調整量として選択できます。調整量の選択は任意であり、結果は変わりません。

idxAdjust = find(labCP.ID==idxInvalid);
labCP.tStop(idxAdjust(1)) = labCP.tStop(idxAdjust(1))-0.5;
labCP.tStart(idxAdjust(2)) = labCP.tStart(idxAdjust(2))-0.5;
labCP(idxAdjust,:)

ans=2×6 table
    ID    tStart    tStop    Censoring    Sex    Lab 
    __    ______    _____    _________    ___    ____

    4         0     49.5         1         0     0.21
    4      49.5       50         0         0      0.2

手順 5. コックス比例ハザードモデルを作成する。

時間に依存しない変数 Sex と時間に依存する変数 Lab があるコックス比例ハザードモデルを当てはめます。

X = [labCP.Sex labCP.Lab];
T = [labCP.tStart labCP.tStop];
b = coxphfit(X,T,'Censoring',labCP.Censoring,'Baseline',0)

コックス比例ハザードモデルを評価する方法の詳細については、打ち切りデータのコックス比例ハザードモデルを参照してください。

参考

coxphfit

共変量が時間に依存するコックス比例ハザード モデル