スプラインをxblcなしで描いてみよう②Coxモデル編

プログラミング

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

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

さて、前回に引き続き、スプラインをxblcコマンド抜きで描いてみよう!というテーマで、今回は生存時間分析によく使われるCox比例ハザードモデルに実装してみたいと思います.

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

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

xblcコマンドはオッズ比やハザード比といった「比」を計算するときに結構トリッキーなことをしてのけていたのですが、以下のような欠点がありました.

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

ハザード比でもオッズ比でも、リファレンスカテゴリを指定する必要があります。そして、指定したリファレンスカテゴリーの値で割り算することで、ハザード比やオッズ比が求められます.

点推定値だけならば話は簡単なのですが、信頼区間があるので、単純に割り算すればよいというわけではありません.

ハザード比やオッズ比の元となる、ハザードやオッズですが、それぞれexp()の中の線形予測子が線形回帰の式になっています.
リファレンスカテゴリのハザードやオッズで割り算する際は、線形予測式の中では引き算になります.

1.ハザード関数

初めに、ハザード比、というもので用いられる、「ハザード」ってなんでしょうか?これはすごく短い時間でイベントが起こる確率のこと、と理解しておきましょう.

もうちょっと正確に言うと、

時間tまで生存していた、という条件の下で、次の瞬間(t+Δt)までの間にイベントが起こる確率

となります.

Cox比例ハザードモデルというのは、ハザードの比を求めることを目的としています.

ロジスティック回帰では確率piのロジット(logit)対数オッズ(log odds) が線形予測式で表現されるのでしたが、ハザードも線形予測式が組み込まれているので、基本は同じです.式にするとこんな感じ.

もう少し詳しく言うと、

ここで少し気を付けたい、というか注目したいのは、β0、つまりy切片がない、というところでしょうか.

それと、線形式の頭にh0(t)というものがくっついていますが、時間tによってウネウネと変化するベースラインハザードになります(ノンパラメトリック).

線形式とノンパラメトリックなベースラインハザードを合わせて使うのでセミパラメトリックモデルと呼ばれるわけです.

共変量z2が1だけ変化したときにハザードがどれだけ変化するかを表したいときは次のようになります.

実際には対数ハザードを計算して最後にexponentして戻してやることになりますが、これはロジスティック回帰と同じですね.

2.スプラインでハザード比を表す

スプラインの場合も、ある曝露因子を連続的に変化させたときに、参照基準におけるハザードと比較して何倍になるのかを見ますので、基本的な考え方としては同じに考えればよいです.

mkspline rcs_idpvar = idpvar, cubic nknot(3)

このように入れると、rcs_idpvar1とrcs_idpvar2という新しい変数が誕生します.

これはそのままでは解読不能なのですが、3次関数の曲線を描くためにいい感じに設定された変数だと思っておきましょう.

stset duration, failure(event==1)
stcox rcs_idpvar1 rcs_idpvar2 cov1 cov2…

ロジスティック回帰の時と同様、残念ながらリファレンスを表現することができません.これをどのように解決するのかを考えてみます.

下の数式はx1をxaとxbの2つの変数に分けたところと考えてください.

リファレンスの位置の rcs_idpvar1, rcs_idpvar2の値はそれぞれxar, xbrとなります.

それをそれぞれxa, xbから引き算すればリファレンスの値を引いた対数ハザード比が出てくるはずです.

引き算した後のスプライン変数をrcs_idpvar1’ rcs_idpvar2’ などとして

係数はそのままにして説明変数 としてそのまま入れて計算させればハザード比がそれぞれ算出されるというわけです.

ここで、信頼区間付きでハザード比を求めるのにやはり”predictnl”を使うわけですが、この時に計算されるのは対数ハザードになります.

なのでここまでの段階でリファレンスの引き算を行ってしまうのが楽です.

3.リファレンスカテゴリのスプライン変数を作成する

ここが今回の記事のメインになります.

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

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

ここで、ある変数Xを横軸にして、3つのノットでリファレンスを10としてRCSを描くプログラムを作ってみましょう.上記の1~3のステップを示しています.

        mkspline X = rcs_index, cubic nknot(3)
* add dummy observation for the reference
	local NEW_OBS = _N + 1
	set obs `NEW_OBS'
	replace X = 10 in `NEW_OBS'
	gen flag=1 if X==0
	sort flag
	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. ハザード比を信頼区間付きで算出する

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

ここで便利な関数コマンド“predictnl”が出てきます.

Cox比例ハザードモデルでの解析を終えた後に、predictnlの関数を使えば対数ハザード比を信頼区間付きで計算してくれます.

predictnl yhat = predict(xb), ci(lb ub) level(95)

  (※lb=lower boundary, ub = upper boundaryのつもり)

あとはyhat, lb, ubをexponentすればOKです.

5.まとめ

ロジスティック回帰のときとはpredictnlに入れるための式が少し異なっていました.

しかしこちらのほうが理解しやすいですね.ロジスティックのときのほうが実は苦労しました.

コメント

  1. stata初心者 より:

    いつも大変参考にさせていただいております。
    質問なのですが、
    mkspline X = rcs_index, cubic nknot(3)
    のrcs_indexは何を指しますでしょうか。
    例えば、血圧(変数:sbp)を横軸にした、予後(変数:event)とのCOX回帰のRCSのグラフを作成したいのですが、
    グラフまでのコマンドを教えて頂けないでしょうか。
    宜しくお願いいたします。

    • 管理者 より:

      ご質問ありがとうございます。
      rcs_indexというのは新しい変数名です。別にどんな名前でもかまわないです。
      Xという変数から、ノットを3つにしてスプライン曲線がかけるように、新たな変数を作ります。
      sbpという変数であれば、スプラインのための変数をsbp_indexとおくと、
      mkspline sbp = sbp_index, cubic nknot(3)
      という風にいれてみてください。あとはここに書かれているようなプログラムを追加していくだけです。

      • stata初心者 より:

        管理者様
        早速のお返事いただき誠にありがとうございます。
        私の勘違いでしたら、大変申し訳無いのですが、
        mkspline sbp_index = sbp, cubic nknot(3)  
        で合ってますでしょうか?
        mkspline sbp = sbp_index, cubic nknot(3)
        ですと、variable sbp_index not found  とでてきます。

        • 管理者 より:

          失礼いたしました。それでできると思います。
          今回のプログラムは煩雑ですので、まずはこちらからお試しいただくといいかと思います。https://statakahiro.com/restricted-cubic-splines%e3%82%92stata%e3%81%a7%e5%ae%9f%e8%a1%8c%e3%81%97%e3%81%a6%e3%81%bf%e3%82%8b

          • stata初心者 より:

            ありがとうございます。そちらも読ませていただき、作成できました。
            多重代入後のHRをスプライン曲線で描きたいと思っております。
            mi xeq: mi estimate を使用するタイミングを教えて頂けないでしょうか。
            よろしくお願いいたします。

  2. 管理者 より:

    ご質問ありがとうございます。
    mksplineのときにまずはmi xeq 0 を頭につける必要があります。
    リファレンスのダミー症例を追加したりのところは不要です。
    モデルを走らせるところは
    mi estimate, post saving(dataset名) : stcox – – – – –
    というのをやったあとに
    mi predictnl b = predict(xb) using dataset名, ci(lb ub)
    という風にしていただく必要があります。ここがComplete-case analysisとの一番大きな違いになります。
    時間を見つけてこれも記事にしていこうかと思います。

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