Poissonでincidence rateを出してスプラインを描く

グラフ

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

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

さて、今回は、非線形の関係にある、連続的な曝露因子に対するIncidence rateを滑らかに表現する方法(スプライン)を示してみたいと思います.

今回の記事を作成するにあたって最初に見た論文がこちら

Teerlink et al. Effect of Ejection Fraction on Clinical Outcomes in Patients Treated With Omecamtiv Mecarbil in GALACTIC-HF. JACC 2021;78:97-108

で、Central illustrationとして描かれていたFigureに興味を持ったからでした.

縦軸をincidence rateとし、横軸は連続値.そして関連は非線形、ということで使い勝手がよさそうです.

しかもよく見れば、このincidence rateを出すのに何やらPoisson regressionを使ったとあるじゃありませんか!謎に満ちていたので、ちょっとやり方をしらべてみました.

1.データセット(オリジナル)の準備

非常に簡単ながら、データセットを用意しておきました.もはやこの手のでたらめなデータセットを作るのは大得意となってしまいましたが、くれぐれも悪用しないでくださいね.間違っても論文投稿などしないようにお願いします.(そんな人いないか…)

さて、データセットです.自己責任でダウンロードお願いします.

これを展開していただくと、771レコードほどのデータが入っていて、case, age, sex, exposure, outcome, year という変数が格納されているので、ご確認ください.

caseとは通し番号で、IDみたいなものです.年齢、性別、曝露因子、アウトカム(0/1)とイベントまでの時間(年)となっています.

2.Poissonでincidence rate ratioを計算する

次に、解析の中心となる、Incidence rateの計算ですが、その前にPoisson回帰でよく求める、Incidence rate ratioを先に求めてみましょう.

import excel using "poisson_rcs.xlsx", firstrow clear 

poisson outcome age sex exposure, exposure(year) irr
ExposureのIRRは0.99と求まった

こんな感じでやることが多いと思います.これをスプラインで表現するとどうなるかが次です.

3.スプラインでIncidence rateを表現する

今度はスプラインを描きたいと思います.最初にノットの数を指定しますが、今回は単純に、3つとしてみますね.

そして、最終的に、incidence rateを100人年あたりで表現したいので、時間(年)を100で割っておきます.

import excel using "poisson_rcs.xlsx", firstrow clear 
mkspline index = exposure, cubic nknots(3) 
gen time = year/100
poisson outcome age sex index* , exposure(time)

このあとは、信頼区間付きでincidence rateを計算するため、Post estimateとして、predictnl コマンドを使用します.

また、Predictするにあたって、covariateすべてに、そのコホートの平均値を入れ替えてください.この作業をしないと一つのグラフになりません.

foreach v of varlist age sex {
    sum `v'
	replace `v' = r(mean)
}	
predictnl ir = predict(ir), ci(lb ub)

このPredictする前の平均への置き換えがなかなか理解されないのですが、しょせんはRCSというのはその集団における代表的な人がExposureだけを変えるとリスクがどんな風に変化するかを見たいだけのものである、ということが根底にあります.

それでも理解できない場合は、そうしないとこうなる、ということだけご理解ください.

わちゃ~

そして無事にPredictができたらtwowayでグラフを描きましょう.incidence rate, 信頼区間、そしてヒストグラムを重ねてみました.

#delimit;
twoway
	(line ir exposure, lc(navy) lp(solid) lw(medthick) sort yaxis(1))
		||
	(line lb exposure, lw(medthick) lc(navy) lp(shortdash) sort yaxis(1))
		||
	(line ub exposure, lw(medthick) lc(navy) lp(shortdash) sort yaxis(1))
		||
	(hist exposure, freq  fcolor(none) lcolor(edkblue) w(10) lw(thin) yaxis(2))
	,
	ytitle("Incidence rate" "(Per 100 Patient-Years)") xtitle("Exposure variable") 
	xscale(range(30 160)) xlabel(50(20)150)
	yscale(range(10 80)) ylabel(0(20)80, ang(0))
	ylabel(0 50 100 150, axis(2)) ysc(range(0 500) axis(2))
	legend(off) graphregion(color(white))
	;
#delimit cr	

そして最後にこんな風にできれば完成です.

まとめ

Poisson regressionからIncidence rateのスプラインを描いてみました.臨床で通常我々が経験するものは、どれも大抵は「珍しいアウトカム」ですから、ほぼどんな場合にも適用可能と考えます.

コメント

  1. […] 以前もxblc抜きでPoisson回帰により算出したincidence rateをスプラインで表現する方法をご紹介しました. […]

タイトルとURLをコピーしました