さまざまなグループのハザード関数と生存時間関数

この例では、さまざまなグループの累積ハザード関数と生存時間関数を推定およびプロットする方法を説明します。

手順 1. 標本データを読み込んで整理する。

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

load('readmissiontimes.mat')

データには、患者の再入院時間と、患者の性別、年齢、体重、喫煙ステータスおよび打ち切りに関する情報が含まれています。このデータは、シミュレーションされたものです。

性別ごとに再入院時間と打ち切りの行列を作成します。

female = [ReadmissionTime(Sex==1),Censored(Sex==1)];
male = [ReadmissionTime(Sex==0),Censored(Sex==0)];

手順 2. 性別ごとに累積分布関数を推定およびプロットする。

女性患者と男性患者について累積分布関数のカプラン・マイヤー推定値をプロットします。

figure()
ecdf(gca,female(:,1),'Censoring',female(:,2));
hold on
[f,x] = ecdf(male(:,1),'Censoring',male(:,2));
stairs(x,f,'--r')
hold off
legend('female','male','Location','SouthEast')

手順 3. 生存時間関数をプロットする。

女性患者と男性患者の生存時間関数を比較します。

figure()
ax1 = gca;
ecdf(ax1,female(:,1),'Censoring',female(:,2),'function','survivor');
hold on
[f,x] = ecdf(male(:,1),'Censoring',male(:,2),'function','survivor');
stairs(x,f,'--r')
legend('female','male')

この図は、男性患者の再入院時間の方が女性患者に比べて短いことを示しています。

手順 4. ワイブル生存時間関数を近似する。

女性患者と男性患者の再入院時間にワイブル分布を近似します。

pd = fitdist(female(:,1),'wbl','Censoring',female(:,2))

pd = 

  WeibullDistribution

  Weibull distribution
    A = 12.5593   [10.749, 14.6745]
    B = 1.99834   [1.56489, 2.55185]

pd2 = fitdist(male(:,1),'wbl','Censoring',male(:,2))

pd2 = 

  WeibullDistribution

  Weibull distribution
    A = 4.63991   [3.91039, 5.50551]
    B = 1.94422   [1.48496, 2.54552]

pd2 = fitdist(male(:,1),'wbl','Censoring',male(:,2))

pd2 = 

  WeibullDistribution

  Weibull distribution
    A = 4.63991   [3.91039, 5.50551]
    B = 1.94422   [1.48496, 2.54552]

推定された生存時間関数に女性患者と男性患者のワイブル生存時間関数をプロットします。

plot(0:1:25,1-cdf('wbl',0:1:25,12.5593,1.99834),'-.')
plot(0:1:25,1-cdf('wbl',0:1:25,4.63991,1.94422),':r')
hold off
legend('Festimated','Mestimated','FWeibull','MWeibull')

ワイブル分布によって、データの適切な近似が得られます。

手順 5. 累積ハザードを推定し、ワイブル累積ハザード関数を近似する。

累積ハザード関数を性別に推定し、ワイブル累積ハザード関数を近似します。

figure()
[f,x] = ecdf(female(:,1),'Censoring',female(:,2),...
'function','cumhazard');
plot(x,f)
hold on
plot(x,cumsum(pdf(pd,x)./(1-cdf(pd,x))),'-.')
[f,x] = ecdf(male(:,1),'Censoring',male(:,2),...
'function','cumhazard');
plot(x,f,'--r')
plot(x,cumsum(pdf(pd2,x)./(1-cdf(pd2,x))),':r')
legend('Festimated','FWeibull','Mestimated','MWeibull',...
'Location','North')

参考

ecdf | coxphfit | ksdensity