2 つのグループの生存時間関数

この例では、2 つのグループに関するデータへのブール型 XII 分布近似を使用して、経験的生存時間関数とパラメトリック生存関数を検出する方法を示します。

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

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

load('lightbulb.mat')

データの 1 列目には 2 種類の電球の寿命 (時間単位) が含まれています。2 列目には、電球のタイプに関する情報が含まれています。0 は蛍光灯電球を、1 は白熱灯電球を示します。3 列目には打ち切り情報が含まれます。1 は打ち切られたデータを示し、0 は正確な故障時間を示します。このデータは、シミュレーションされたものです。

電球タイプごとに 1 つの変数を作成し、打ち切り情報も含めます。

fluo = [lightbulb(lightbulb(:,2)==0,1),...
			lightbulb(lightbulb(:,2)==0,3)];
insc = [lightbulb(lightbulb(:,2)==1,1),...
			lightbulb(lightbulb(:,2)==1,3)];

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

2 種類の電球の推定生存時間関数をプロットします。

figure()
[f,x,flow,fup] = ecdf(fluo(:,1),'censoring',fluo(:,2),...
				'function','survivor');
ax1 = stairs(x,f);
hold on
stairs(x,flow,':')
stairs(x,fup,':')
[f,x,flow,fup] = ecdf(insc(:,1),'censoring',insc(:,2),...
				'function','survivor');
ax2 = stairs(x,f,'color','r');
stairs(x,flow,':r')
stairs(x,fup,':r')
legend([ax1,ax2],{'Fluorescent','Incandescent'})
xlabel('Lifetime (hours)')
ylabel('Survival probability')

Figure contains an axes object. The axes object with xlabel Lifetime (hours), ylabel Survival probability contains 6 objects of type stair. These objects represent Fluorescent, Incandescent.

白熱灯電球の生存確率が、蛍光灯電球の生存確率よりも大幅に低いことがわかります。

手順 3. ブール型 XII 分布を近似する。

ブール分布を蛍光灯タイプの電球と白熱灯タイプの電球の寿命データに近似します。

pd = fitdist(fluo(:,1),'burr','Censoring',fluo(:,2))

pd = 
  BurrDistribution

  Burr distribution
    alpha = 29143.4   [0.903922, 9.39617e+08]
        c = 3.44582   [2.13013, 5.57417]
        k = 33.7039   [8.10737e-14, 1.40114e+16]

pd2 = fitdist(insc(:,1),'burr','Censoring',insc(:,2))

pd2 = 
  BurrDistribution

  Burr distribution
    alpha = 2650.76   [430.773, 16311.4]
        c = 3.41898   [2.16794, 5.39197]
        k =  4.5891   [0.0307809, 684.185]

ブール型 XII 生存時間関数を重ね合わせます。

ax3 = plot(0:500:15000,1-cdf('burr',0:500:15000,29143.5,...
			3.44582,33.704),'m');
ax4 = plot(0:500:5000,1-cdf('burr',0:500:5000,2650.76,...
			3.41898,4.5891),'g');
legend([ax1;ax2;ax3;ax4],'Festimate','Iestimate','FBurr','IBurr')

Figure contains an axes object. The axes object with xlabel Lifetime (hours), ylabel Survival probability contains 8 objects of type stair, line. These objects represent Festimate, Iestimate, FBurr, IBurr.

ブール分布は、この例の電球の寿命に良好な近似を提供します。

手順 4. コックス比例ハザードモデルを当てはめる。

電球のタイプが説明変数であるコックス比例ハザード回帰を近似します。

[b,logl,H,stats] = coxphfit(lightbulb(:,2),lightbulb(:,1),...
'Censoring',lightbulb(:,3));
stats

stats = struct with fields:
                    covb: 1.0757
                    beta: 4.7262
                      se: 1.0372
                       z: 4.5568
                       p: 5.1936e-06
                   csres: [100x1 double]
                  devres: [100x1 double]
                 martres: [100x1 double]
                  schres: [100x1 double]
                 sschres: [100x1 double]
                  scores: [100x1 double]
                 sscores: [100x1 double]
    LikelihoodRatioTestP: 0

$p$ 値 p は、電球のタイプが統計的に有意であることを示しています。ハザード率の推定値は $e x p$ ( $b$ ) = 112.8646 です。つまり、白熱灯電球のハザードは、蛍光灯電球のハザードの 112.86 倍であることを意味します。