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

カテゴリ変数について割合の信頼区間を計算するには、通常、二項分布を仮定してClopper and Pearsonの方法を使う（『Rによる保健医療データ解析演習』の第9章、9.3に書いた通り）。しかし、カテゴリが3つ以上ある場合、各カテゴリの割合は独立ではない（もちろん2つの場合も片方のカテゴリの割合が決まれば、残りも自動的に決まってしまうわけだが、1つのカテゴリが生起する確率、と考えれば二項分布を普通に適用して差し支えない）ので、別のアプローチをした方が良いと思われた。

自分の本にも書かなかったけれども、カテゴリ３つ以上の場合の、各カテゴリの割合の信頼区間の推定についても、既に多くの研究があり、Rでも実装しているパッケージがある。DescToolsパッケージのMutinomCI()関数である（もっとも、DescToolsパッケージのvignettesの１つの13.1の解説にあるように、カテゴリ３つ以上でも二項分布を当てはめて信頼区間を求める、BinomCI()という関数も用意されているので、必ずしも各カテゴリについて二項分布で信頼区間を推定して悪いわけではないと思われる）。

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.99.60を文献として引用すれば良く、カテゴリ３つ以上の信頼区間の推定方法としては、Sison and Glaz (1995)の方法を使ったと書いて、Sison CP and Glaz J (1995) Simultaneous confidence intervals and sample size determination for multinomial proportions. Journal of the American Statistical Association, 90:366-369.を引用すれば良い。

使い方は簡単で、例えば、『Rによる保健医療データ解析演習』の9.5.1で、1日当たりの交通事故件数の分布がポアソン分布に従うかどうかをカイ二乗適合度検定で調べる方法を説明したデータを使えば、まず、install.packages("DescTools", dep=TRUE)でDescToolsパッケージをインストールしておいて、それをlibrary(DescTools)で呼び出したら、後はいきなりMultinomCI(c(79, 61, 13, 1, 1), method="sisonglaz", conf.level=0.95, sides="two.sided")と打てば良い。95％両側信頼区間をSison and Glaz (1995)の方法で求めるというデフォルトで良いならばこれらのオプション指定は不要なので、たんにMultinomCI(c(79, 61, 13, 1, 1))で良い。結果として返ってくるオブジェクトは、列名付きの行列なので、例えばresというオブジェクトに付値しておけば、sapply(as.data.frame(res), function(x) sprintf("%4.3f", x))によって、小数点以下3桁だけ表示させることもできる。

コメント

コメントを投稿

このブログの人気の投稿

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=...
続きを読む

オッズ比よりもリスク比を推定すべきという話

 院生から相談があったわけではないが、何回かのメモから採録する。    【第1807回】 国際保健医療学会西日本地方会など（2025年3月1-2日） いくつかの発表について三重大の谷村先生が指摘していたのは、 NEJMのStat Guideline にも書かれている prevalenceが高いときにはORは過大評価になるのでロジスティック回帰でなく修正ポアソン回帰を使ってCPRを出すべき （リンク先のプレゼン資料が大変わかりやすい）という話。参考に "Odds Ratio or Prevalence Ratio? An Overview of Reported Statistical Methods and Appropriateness of Interpretations in Cross-sectional Studies with Dichotomous Outcomes in Veterinary Medicine" (Frontiers in Veterinary Science, 2017)、 "Commentary: Calculating risk and prevalence ratios and differences in R: developing intuition with a hands-on tutorial and code" (Annals of Epidemiology, 2023)、 "Incorrect inference in prevalence trend analysis due to misuse of the odds ratio" (Annals of Epidemiology, 2016)、 "Overestimation of Relative Risk and Prevalence Ratio: Misuse of Logistic Modeling" (Diagnostics, 2022)、辺りは読んでおくべきか。ちなみに修正ポアソン回帰は、Rならsandwichパッケージかbootパッケージかrqlmパッケージを使うのが便利そう。   https://minato.sip21c.org/im3r/2025030...
続きを読む

対応のある多群間のノンパラメトリックな比較のための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弾）===...
続きを読む