スプラインをxblcなしで描いてみよう③異なるサブグループに分けた曲線を重ねる

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

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

さて、前回に引き続き、スプラインをxblcコマンド抜きで描いてみよう！というテーマで、今回は生存時間分析によく使われるCox比例ハザードモデルを基に、複数のサブグループのグラフを重ねるグラフを描いてみます．

ロジスティック回帰のスプラインはこちら

Poisson回帰によるincidence rateのスプラインはこちら

Cox回帰のスプラインはこちら

xblcコマンドは初学者にとっては非常に使いやすくて便利なコマンドなのですが、以下のような欠点がありました．

異なる複数群で別々のRCSを引くときに、リファレンスを片側の群で設定することができない
多重補完後のデータで行うことができない
やたらと時間がかかる

サブグループ解析などでよくみられるのですが、異なる複数群で、共通したリファレンスだけど別々の曲線を重ねる、という方法が必要なことがあります．

例として挙げるのは、こちらの論文です．

Grams M, et al. A Meta-analysis of the Association of Estimated GFR, Albuminuria, Age, Race, and Sex With Acute Kidney Injury. Am J Kidney Dis. 2015 Oct; 66(4): 591–601.

これは、メタ解析ではありますが、例えばFigure 1では、eGFR 80, アルブミン尿なしの人をリファレンスとして、Moderately-increased UACR, Severely-increased UACRそれぞれの群で、eGFR全域と比較したハザード比を出していますね．

xblcではこのようなことはできないのです．ではどうすればいいのか？

xblcを使わないスプライン曲線の作成に慣れることです．

1．手順①
4．手順②
5．まとめ

1．手順①

別々のサブグループの結果を重ねるので、単純にグループ変数とifを使って別々に解析をすればよいのか、というとそうではありません．リファレンスを共通にするということはそのグループを含んだ回帰式を作らねばならないのです．

そこで重要なのは、グループ変数とスプライン変数との交互作用項を入れることです．

手順としては次のように行います．

テーブルの一番下に新しいダミーレコードを作成する
ダミーレコードにおいて、RCSの横軸の変数Xがリファレンスとなる値を、その変数Xの中に代入する
変数Xがダミーレコードと同じ症例から、スプライン変数を取得する
新たに生成したスプライン変数とグループ変数の交互作用項を組み込んでCox比例ハザードモデルを実行
グループ変数の各群における点推定値と信頼区間をpredictnlで算出する
グラフの作成

ここで、ある変数Xを横軸にして、3つのノットでリファレンスを10としてRCSを描くプログラムを作ってみましょう．上記の1～3のステップを示しています．（ここまでは以前の記事と全く同じです）

        mkspline rcs_index = X, cubic nknot(3)
* add dummy observation for the reference
	gen flag=1 if X== 10
　　　　sort flag

　　　　local NEW_OBS = _N + 1
	set obs `NEW_OBS'
	replace X = 10 in `NEW_OBS'	
	local rcs_ref1 = rcs_index1[1]
	local rcs_ref2 = rcs_index2[1]

	replace rcs_index1 = `rcs_ref1' in `NEW_OBS' 
	replace rcs_index2 = `rcs_ref2' in `NEW_OBS'
	
forvalue k = 1/2 {
	local RCS_REF`k' = rcs_index`k'[_N]
	disp `RCS_REF`k''
}	
* comparing to the reference
forvalue k = 1/2{
		local RCS_REF`k' = rcs_index`k'[_N]
		replace rcs_index`k' = rcs_index`k' - `RCS_REF`k''
	}

4．手順②

次に、ハザード比を信頼区間付きで算出する手順を説明していきます．

以前の記事では、Cox回帰の部分は端折ってしまいましたが、今回はここからが肝ですので、端折らずに示します．グループ変数をgroupとします．（1.アルブミン尿なし、2.中等度、3.重度、と解釈してもらえればいいと思います．）

foreach var of varlist covariates {
  sum `var', meanonly
  replace `var' = `var' - r(mean)
  }  …⓪
stcox c.rcs_index*##i.group covariates
　…①
test c.rcs_index1#2.group c.rcs_index2#2.group c.rcs_index1#3.group c.rcs_index2#3.group
　…②
foreach v of varlist covariates {
  sum `v’
　replace `v’ = r(mean)
　｝
　…③
forvalues n = 1/3 {　replace group = `n’
  predictnl yhat`n’ = predict(xb), ci(lb`n’ ub`n’)  }
　…④

1行ずつ解説してみたいと思います．

⓪　中心化　(2022/5/28追記)

foreach var of varlist covariates {
  sum `var', meanonly
  replace `var' = `var' - r(mean)
  }

イタリックで記載されたcovariatesは、回帰式にいれるための共変量です．いくつ入れてもいいのですが、交絡因子と思われるものを入れていきます．

こちらのステップは回帰分析の結果自体には実は何ら影響しませんが、predictするときに影響があるようで、結論から言えば毎回このステップを踏んでおいた方が安全です．

中心化の前後で、ハザード比の点推定値、信頼区間ともに全く同じ値になるので影響なさそうなのですが、最後のステップですべてを平均に置き換える段階（ステップ③）で、中心化をしたかしていないかによって結果が変わります．

もしこの段階で中心化をしてない場合には③はr(mean)ではなく、0に置き換えるようにしてください．

①　回帰式の実施

stcox c.rcs_index*##i.group covariates

まず最初のCoxを実行するところで、スプライン変数であるrcs_index1, rcs_index2という連続変数と、グループ変数との間の交互作用項を入れているのがわかりますか？

グループ変数はカテゴリ変数ですので、必ず「i.」を先頭に入れましょう．参照カテゴリは一番小さいところに置かれますが、ib2, ib3などとやれば参照カテゴリは自由に変えることができます．

②　交互作用の検定

test c.rcs_index1#2.group c.rcs_index2#2.group c.rcs_index1#3.group c.rcs_index2#3.group

次の行は、この交互作用項が有意かどうかの検定をしています．グループ変数リファレンスがグループ1ですので、グループの数がk個あるとしたら、以下のようになります．

group2×スプライン変数1　group2×スプライン変数2　… groupk×スプライン変数1　groupk×スプライン変数2

といった具合に並べていくようにしてください．スプラインのノットの数によってはスプライン変数の数が増えますので、よりたくさん書かないといけませんね.

③　共変量を揃える

foreach v of varlist covariates {
  sum `v’
　replace `v’ = r(mean)
　｝

さて、その次ですが、ここは共変量をすべて同じ値に揃えます．そうしないとグラフがガタガタになるという話は以前した通りです．

繰り返しになりますが、もしステップ⓪で中心化をしてない場合にはここはr(mean)ではなく、0に置き換えるようにしてください．

④　グループごとに予測値を求める

forvalues n = 1/3 {  replace group = `n’  predictnl yhat`n’ = predict(xb), ci(lb`n’ ub`n’)  }

ここが一番の肝になるのですが、

一旦すべての人をそれぞれのグループの人にすり替えてしまいます．

それらの人たちはすでに共変量が共通したものになっていますので、このようなことをしても問題ないわけです．なんだか不思議な感じがするかもしれませんが．

それからこの値をexp()で指数化し、HRと信頼区間を得ます．それを基にして、twowayのグラフを作成して重ね合わせましょう．

5．まとめ

サブグループごとのスプライン曲線を、共通したリファレンスをもって重ね合わせるための方法をご紹介しました．

結構この方法は需要があると思いますので、多くの方に役立てばよいと思っております．

1．手順①

4． 手順②

⓪ 中心化 (2022/5/28追記)

① 回帰式の実施

② 交互作用の検定

③ 共変量を揃える

④ グループごとに予測値を求める