比例ハザード性からの逸脱への対応

このブログでは、統計解析ソフトStataのプログラミングのTipsや便利コマンドを紹介しています．

Facebook groupでは、ちょっとした疑問や気づいたことなどを共有して貰うフォーラムになっています． ブログと合わせて個人の学習に役立てて貰えれば幸いです．

Cox比例ハザードモデルは、生存時間分析を実施する上で最もオーソドックスな統計手法ですが、これを使うための最低条件として、

「比例ハザード性」

というのがあります.

これはどういうことかというと、カプランマイヤー曲線を思い浮かべていただくとわかりやすいのですが、ある状態Aと状態Bにおいて、無イベント生存率が平行して下がっていくような状況においてのみ使えますよ、ということです．

比例ハザード性からの逸脱を客観的に示す方法としては、Graphicalに示す方法と、統計学的な手法とがあります.

前者はlog-log plotと呼ばれる方法が良くとられますが、Schoenfeld残差を利用した方法もあります．また同様にSchoenfeld残差を利用した統計学的な方法もあります．これらは前回の記事に掲載しております．

そして、これらの方法で比例ハザード性がよさそうです、ということになれば晴れてCoxを使えるわけですが、時々問題あり！となることがあります．

そこで、今回は、そういった状況においてどのように対処するかをまとめてみました．

1．Piecewiseにハザード比を求める

この手法はフォロー期間を細かく区切ってその中でのハザード比を出すというものです．独立変数の影響の経時的な変化を出すのに有効です.

区切れ目と水準の数を決めるのが最初のステップになります．

透析患者と各種固形癌の予後を比べるという研究にちょうどよい例がありましたのでご紹介します．カプランマイヤーからして明らかに比例ハザード性からの逸脱が見られそうです．実際、イベントの起こるメカニズムが異なるとき必ずしも平行にならないということは以前から知られていました．

そんな前提なので区間ごとにHRを出すというこのやり方をとって、1年ごとのHRの推移を図にしていたのでした．

こちらの研究も同様に区間ごとにHRをみていますが、この研究では、あるマーカーのアウトカムへの関連が時間の経過とともに薄まってしまう、という結果になりました．そういった関係を可視化する上でも有効な手法になります．

Stataのコマンドではstsplitというのがあります．次のデータセットでまず比例ハザード性の確認をします．

use http://www.biostatepi.org/data/hers, clear
stset pafu, f(pa) id(id)
   /* idの指定を忘れずに */
stcox group
estat phtest, plot(group)
   /* 比例ハザード性逸脱を確認 */  

1年ごとに区切ってみましょう．そこでstsplitの出番になります．

stsplit dur, at(0(1)5) 
/* 期間を1年ごとに区切ってみます */	
forvalues n = 0/5 {
	stcox group if dur==`n'
	}
	/* あるいはループを使って一気に答えをだすと */
qui {
forvalues n = 0/5 {
	stcox group if dur == `n'
	local coef_`n' = _b[group]
	local low_`n'  = _b[group] - invnormal(0.975)*(_se[group])
	local high_`n' = _b[group] + invnormal(0.975)*(_se[group])
noi disp %3.2f exp(`coef_`n''), %3.2f exp(`low_`n''), %3.2f exp(`high_`n'')
	}
}
	/* この結果をもって作図すると */
preserve
clear
input hr low high time
1.52 1.01 2.29 0
0.98 0.66 1.46 1 
0.85 0.54 1.33 2 
0.64 0.39 1.06 3
1.03 0.50 2.14 4
1.00 1.00 1.00 5
end

#delimit;
	twoway rcap high low time, lcolor(maroon)
	|| 
	connected hr time, sort msymbol(Dh) lcolor(maroon) mcolor(maroon)
	,
	xlabel(0(1)5) legend(off) 
	yline(1, lcolor(black)) yscale(range(0.5 2)) 
	ylabel(0.5(0.5)2, angle(0)) ytitle("Hazard ratios") 
	xtitle("Observation period") 
	;
#delimit cr		
restore  

といった具合に書くと次のようなグラフになります．

まあちょっと荒々しいですが、最初の1年間はハザードが高いけど、徐々に下がっていくというようなことがわかりますね．

2．時間×曝露因子の交互作用項を入れる

これはメインのアウトカムとして示す以外に、時間依存変数かどうかを検証するのにも用いられる方法です（昨日の記事）．

stcox varlist, tvc(var) texp(ln(_t))

といった形で示すことが一般的のようです．なぜここで時間の自然対数を入れるのか？と疑問に思われた方もいると思いますが、これはフォローアップ時間の分布がどうしても右に細長く尾を引いてしまうためです．

UCLAのサイトでも紹介されています．こちらの動画でも同様に時間の自然対数を入れています．

同じデータセットを使って実際にtvcを走らせてみると、以下のような結果が出力されます.

tvcのところのgroupをみると確かにP<0.05となっていますので、時間依存性であると結論付けることができるようです．

3．スプラインで滑らかに示す

ぶっちゃけていってしまえば、前二つは前座で、今日のメインはこちらです．スプラインにしてしまうという方法です．これは１の方法の拡張版とみなすこともできるかもしれません．なめらかなので、信頼区間もある程度限定されますし、信頼区間があるからこそ細かい時間で検証できるので、重宝しそうです．

コマンドは“stphcoxrcs”というそのものズバリな名称です．

使用上のポイントは3つです．

1. stsetの段階で、id()で個人を指定しておく必要がある
2. Coxモデルを走らせた直後にコマンドを実行するだけ
3. オプションにtwoway optionがあるのでグラフを重ねることができる

use http://www.biostatepi.org/data/hers, clear
stset pafu, f(pa) id(id)
   /* idの指定を忘れずに */
stcox group
estat phtest
   /* 比例ハザード性逸脱を確認 */
stphcoxrcs group, splitevery(.2) nknots(3) range(.25 6) noyref

これで出てくるのが下の図です．

大体時間0.5～2くらいの間はHRが相対的に高い時期であることがわかります.

これは１のところで飛び飛びの値で作ったグラフよりもなめらかですし、信頼区間が狭い分美しいですよね．

こういった情報を可視化できる強力なツールといえると思います．

コマンドの詳細はhelp stphcoxrcsで調べてみてください．

4．まとめ

比例ハザード性からの逸脱が見られた時の対処方法についてまとめてみました．

ちなみに以前一緒に働いていたPennsylvania大学の生物統計の専門家は、比例ハザードについては図と統計の両方を一応確認するけど、明らかにおかしい、というもの以外は特別な対処をしないといっていました．サンプル数が増えれば有意になるし、無理に調整しようとして複雑になりすぎて結果的に解釈困難、説明困難な複雑怪奇なモデルになってしまうのは避けたほうが良い、ということでした．何とも現実主義的ですね．極端な原理主義に走る必要はないよ、ということかなと思います．諸説あるでしょうが．

また、今一緒にパートナーを組んで仕事させていただいている統計家の先生も似たようなご意見でした．臨床疫学は臨床家が読者ですから、読み手フレンドリーなほうがいいでしょうね．