SASではcompeting riskの解析方法が今まで正式にサポートされてこなかったため、Gray (1988)の方法を実装した%CIFマクロを使うしかありませんでした。 ちなみに、Gray (1988)はKaplan–Meier推定量ベースのCIF推定量と、それを利用したclass K統計量を提案しており、ノンパラメトリックなcompeting riskの解析方法として広く利用されてきました。
最近になって SAS/STAT 13.1で機能追加が行われ、Procedureだけで基本的なcompeting riskの解析を全て行う事ができるようになりました。 SAS/STAT 13.1以降では、PHREG Procedureのmodel statementでeventcodeオプションを指定することで、Fine and Gray (1999)のモデルが利用できるようになりました。 Fine and Gray (1999)のモデルも、Coxモデルベースのcompeting riskの解析方法として広く利用されていると思います。 この機能追加にともなって、baseline statementで、Fine and Gray (1999)のbaseline cumulative subdistribution hazardの推定量が計算できるようになりました。 Baseline cumulative subdistribution hazardは、Breslow推定量ベースのCIF推定量であり、共変量がない場合にGray (1988)のCIF推定量とほぼ同じ推測が可能ではないかと思います。 それぞれKaplan–Meier推定量ベースとBreslow推定量ベースの方法であるため、数値的には若干結果が異なりますが、両者は漸近的に同等であるため、SASでも簡単にCIFの計算ができるようになったと言えるでしょう。
cif.sas (サンプルプログラム)
SAS/STAT 14.1でさらなる機能追加が行われ、LIFETEST Procedureでノンパラメトリックなcompeting riskの解析(%CIFマクロに相当)が実施できるようになりました。 Gray (1988)の検定やノンパラメトリックなCIF推定量が計算できます。
cif_lifetest.sas (サンプルプログラム; LIFETEST Procedure; SAS/STAT 14.1以降で利用可能)
上記のサンプルプログラムを用いてPHREG Procedureで解析を行った結果と、サンプルプログラムには含めていませんが%CIFマクロで解析を行った結果を以下に示します。 CIFの群間差の検定結果も、CIFの推定結果も大きな差はありませんでした。
PHREG Procedureの検定結果です (STAT 13.1では検定はWald型のみ出力可能です)。
Type 3 Tests Wald Effect DF Chi-Square Pr > ChiSq group 1 3.1487 0.0760
%CIFマクロの検定結果です (Gray's testの結果; class K統計量なので、Coxモデルのスコア型の検定に近い)。
Gray's Test for Equality of Cumulative Incidence Functions Chi- Pr > Square DF Chi-Square 3.10766 1 0.0779
LIFETEST Procedureの検定結果です (SAS/STAT 14.1以降で利用可能)。
Gray's Test for Equality of Cumulative Incidence Functions Pr > Chi-Square DF Chi-Square 3.1077 1 0.0779
PHREG Procedureのbaseline cumulative subdistribution hazardの推定結果をグラフにしたものです。標準誤差の推定方法はシミュレーションベースで実装されており、サンプルプログラムでは25000回反復して求めています。
LIFETEST ProcedureのCIF推定結果をグラフにしたものです。
%CIFとLIFETESTは完全に一致し、PHREGとその他は多少違いがありますがほとんど同じグラフが描けました。
proc phreg data = example; class group / param = ref ref = first; model time*event(0) = group / eventcode = 1; run;
proc phreg data = example; model time*event(0) = / eventcode = 1; strata group; baseline out = out1 cif = _all_ / normalsample = 25000 cltype = loglog; run;
proc lifetest data = example outcif = out1; time time*event(0) / eventcode = 1; strata group; run;