区間打ち切りデータのコックス回帰

    • 駒場時代からの友人であるA君からRで区間打ち切りのある生存時間データをコックス回帰する方法について相談を受けたので、少し調べてみた。

      彼は区間打ち切りデータ（interval censored data）を解析する（icenRegパッケージ）[R]（ねこすたっと）を参照し、icenRegパッケージを使おうとしてエラーに見舞われたとのこと。

      とりあえず手元で動くかどうかを試してみたところ（ただしtidyverseは嫌いなので、tidyverseを使わない書き方[d <- IR_diabetes; d$left <- ifelse(d$left < quantile(d$left, 0.1), 0, d$left); d$right <- ifelse(d$right > quantile(d$right, 0.9), Inf, d$right);]でデータフレームdを作成したが）、core i9マシンだというのにかなり時間は掛かったものの、無事に動いたことから、たぶんパッケージインストール時に依存パッケージを自動的にインストールするdep=TRUEオプションを付けなかったのではないかと思われたので、そのようにアドバイスした。

      区間打ち切りデータへのコックス回帰は、Rに標準実装されているsurvivalパッケージのcoxph()ではできないのだが（注：カプラン＝マイヤ法で生存時間の中央値を推定するsurvfit()や、パラメトリックモデルを当てはめるsurvreg()は区間打ち切りデータも使える）、cranを探ると、何もicenRegに限らず、他にもいくつかパッケージがある。ICsurvパッケージやICcalibパッケージ（論文も出ている）、intccrパッケージが見つかった。それぞれ癖があって既存データに簡単に使うことは難しかったが、もしパッケージ依存以外の理由でicenRegパッケージが動作しないなら、これらを試してみるのも良いかもしれない。

      あと、SmoothHazardパッケージも使えるかと思ったが、これは加速モデルのようだった。また、元のコックス回帰と違いがあるのかもしれないが、区間打ち切りデータのためのCox-Aalenモデルも、既にcranからは除去されてしまっているが、coxintervalパッケージにcoxaalen()という関数として実装されていた。

      昔タシンボコ研究で出生分析をしたとき（人口学研究に2000年に研究ノートとして載った）、出生間隔のデータに区間打ち切りが入っていて、当時は東大にいたのでSASで分析した（第1子と第2子の出産間隔の分析だったので、産後無月経期間を考慮したモデルについてはCでコードを書いた）ことがあったので、あれをRでやり直してみるかな。というわけで、BI_TF.Rを書いてみた。

      ちなみに、RのSurv(start, end, (event, )type="right/left/counting/interval/interval2")で区間打ち切りを示すためにはイベントフラグ（event）を3にするが（イベントフラグが3でないレコードについてはendは無視される）、区間の与え方にはtype="interval"とtype="interval2"の2種類がある。前者が(start, end]であって、event=を引数として与えねばならないのに対して、後者では、左側打ち切り（event=2）が(-Inf, end)、右側打ち切り（event=0）が(start, Inf)、正確なイベント発生時（event=1）は(start, end)（注：この時はstart=endのはずだから、[start, end]と書くべきかもしれないが、helpには(t, t)と書かれている）、区間打ち切り（event=3）が(start, end)を示し、type="interval2"の方が役に立つことが多いとわかってきたとhelpに書かれている。なお、無限を示すにはInfでもNAでも良いとのこと。重要なポイントは、type="interval2"の場合、startとendの値から自動的にeventは判定されるので、引数としてeventを与えてはいけない（エラーになる）ということである。

      A君から再びメールがあり、今度はflexsurvパッケージのflexsurvspline()関数を試しているそうだ。Rはできることは多いのだが、パッケージの種類が多すぎ、それぞれ癖があるのが初心者に優しくないところだと思う。できるだけ追加パッケージを使わなくてもできる解析にしておく方が安全なのだが、やりたいことが複雑だとそうもいかないのが辛いところ。

      メールの返事に書いた参考情報をここにも採録しておこう。

      Rの関数のヘルプは、例えばicenRegパッケージに入っているic_sp()関数のヘルプを見たいなら、library(icenReg); ?ic_spで良い。もっと見やすく表示されるサイトがいくつかあって、例えば、ic_sp()を見るには、r-universeや、r-packages.ioの該当ページ（検索もできる）が見やすい。ついでにr-universeで見つけたicsp2パッケージのic_sp2()関数は区間打ち切りデータに対してコックス回帰をするのに良いかもしれない。統計学的に高度なことをやっているパッケージの多くは数式の丁寧な説明が論文になっているので、それを読んでおく方が良く、icenRegも論文になっている。先人のQ&AがStack ExchangeやStack Overflow（リンク先はともにic_spを検索した結果）に蓄積されているので、それが参考になることも多い。どうしても困った時はr-helpというメーリングリストがあるので、そこで再現可能なコード例を付けて尋ねると、大抵の場合、誰か答えてくれる人がいるのだが、helpを見ればわかるようなことを尋ねると叱られるので要注意（自分が書いた2冊目のR本の1.1.1 フリーソフト利用上の心得に書いた）。 


    • コメント

    コメントを投稿

    このブログの人気の投稿

    対応のある多群間のノンパラメトリックな比較のためのFriedmanの検定後の対比較について

    Friedmanの検定後のpost hocな多重比較の結果がEZRとSPSSとjamoviで食い違うのは何故かという問い合わせメールをもらったので調べてみた。以下メールの返事を貼り付けて若干改変。 （返信第1弾）============================= jamoviは、 https://bookdown.org/sbtseiji/jamovi_complete_guide/sec-ANOVA-friedman.html に書かれている通り、Friedmanの後の対比較ではbonferroni補正なしのDurbin-Conover法を使っています。  このDurbin-Conover法というのは検定の多重性の調整をしていません。内部的にはPMCMRplusというパッケージに含まれている、frdAllPairsConoverTest()という関数を使っているようで、 この関数自体にはオプションでp.adj="bon"を指定できる のですが、jamoviではそこをいじれない実装になっています。 もっとも、Bonferroniの補正という手法は、対比較で出てきたp値に比較回数を掛けるだけなので、5群間の2群ずつの総当たり比較なら10回の比較がされているので、得られたp値を10倍すれば良いわけです。 （データを使った検証部分は省略） EZRが何故合わないかですが、EZRは対比較にWilcoxonの 符号順位検定を使っています（とメッセージにでてきます）。EZRでは正規近似でp値を出しているので、同順位の値が含まれていると、 wilcox.test.default(pairwise.response1, pairwise.response2, alternative = "two.sided", で警告がありました: タイがあるため、正確な p 値を計算することができません という警告が出ます。 ここで表示されているp値は、wilcox.test(..., paired=TRUE)で得られるp値の10倍になっています。同順位があっても正確なp値が欲しい場合は、exactRankTestsパッケージのwilcox.exact()関数を使えばOKです（もちろん個々のp値は10倍する必要があります）。 （返信第2弾）===...
    続きを読む

    カテゴリ3つ以上の割合の信頼区間

    カテゴリ変数について割合の信頼区間を計算するには、通常、二項分布を仮定してClopper and Pearsonの方法を使う（ 『Rによる保健医療データ解析演習』 の第9章、9.3に書いた通り）。しかし、カテゴリが3つ以上ある場合、各カテゴリの割合は独立ではない（もちろん2つの場合も片方のカテゴリの割合が決まれば、残りも自動的に決まってしまうわけだが、1つのカテゴリが生起する確率、と考えれば二項分布を普通に適用して差し支えない）ので、別のアプローチをした方が良いと思われた。 自分の本にも書かなかったけれども、カテゴリ３つ以上の場合の、各カテゴリの割合の信頼区間の推定についても、既に多くの研究があり、Rでも実装しているパッケージがある。 DescToolsパッケージのMutinomCI()関数 である（もっとも、 DescToolsパッケージのvignettesの１つ の13.1の解説にあるように、カテゴリ３つ以上でも二項分布を当てはめて信頼区間を求める、BinomCI()という関数も用意されているので、必ずしも各カテゴリについて二項分布で信頼区間を推定して悪いわけではないと思われる→正しくは、個々のカテゴリの割合についての95％信頼区間を求めるのか、すべてのカテゴリの割合について95％同時信頼区間を求めるのか、という違いで、どちらも間違いではない。黒木さん、ご指摘ありがとうございました）。 MultinomCI()関数にはさまざまな信頼区間推定方法が実装されていて、古くからあり他のパッケージでも実装されているものもあるらしいGoodman (1965)の方法も、method="goodman"というオプションを指定すれば使えるが、デフォルトがSison and Glaz (1995)の方法（明示的にオプション指定するにはmethod="sisonglaz"だが、省略できる）なので、特別な理由がなければそれで計算するのが良いと思う。DescToolsパッケージを引用するには、Signorell A (2025). _DescTools: Tools for Descriptive Statistics_.  doi:10.32614/CRAN.package.DescTools , R package version 0...
    続きを読む

    CochranのQ検定の話

    以前投稿した、 https://stat-consul.blogspot.com/2024/03/friedman.html とも少し関連しているが、 対応のある3群以上の2値データの比較の話で相談に乗ったのでメモ。  jamoviを使っている院生からの相談に関連して。2時点間でMcNemarの検定をするような、0/1のカテゴリ変数の割合の変化がランダムではない一方向への変化なのかを、3時点以上で検定する方法としてはCochranのQ検定を使うのが普通だが、 jamoviでは特別なメニューが用意されていないので、Friedmanの検定をすればCochranのQ検定と同じp値が得られる 、という記述を見つけたので、確かめてみる。Rのパッケージとしては、RVAideMemoireパッケージに入っているcochran.qtest()関数が便利なので（DescToolsパッケージのCochranQTest()関数もほぼ同様である）、そのexample()の例で試してみると、 response <- c(0,1,1,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,1,0,0,1) fact <- gl(3,1,30,labels=LETTERS[1:3]) block <- gl(10,3,labels=letters[1:10]) library(RVAideMemoire) cochran.qtest(response ~ fact | block) とすれば、 Q = 10.8889, df = 2, p-value = 0.00432 というCochranのQ検定の結果、A、B、Cそれぞれにおける1の割合、AとB、BとC、CとAという対比較をfdr法で検定の多重性を調整した結果のp値が示される。反復測定デザインとしては、blockが個人、factが測定日に当たると考えれば良い。このデータをデータフレームddにしてタブ区切りテキスト形式でファイル出力するには、 dd <- data.frame(response=response, fact=fact, block=block) write.table(dd, "cochran.txt", row.names=FALSE, sep=...
    続きを読む