副作用のためのベータ分布

副作用予測のためのベータ分布 医学のための統計学

• 管理をkindleにやってもらってみた(上記表紙とこちら)。値段を最低価格でつけないといけないのは、epub形式管理の代償ということで…
• epub化はこちらやこれを使う

Rstudio/knitr/pandocでepubを作る

• 表紙の作り方はこちら
• まずはRコマンドのみ

N1 <- 1000
S1 <- 4
N2 <- 10
S2 <- 0
# K/Nの値の想定値をいろいろとる
p <- seq(from=0,to=1,length=10000)
# ベータ分布に準えて、(N,S)の値をあてはめる
p1 <- dbeta(p,S1+1,N1-S1+1)
p2 <- dbeta(p,S2+1,N2-S2+1)
matplot(p[1:50],cbind(p1,p2)[1:50,],type="l")
# それぞれの次の利用時の事件の発生確率
(S1+1)/(N1+2)
(S2+1)/(N2+2)

# riskの下駄
# K/Nの値は0-1でどれかわからない、ということはないでしょう、かなり0に近いでしょう、というときは下駄を履かせます
a <- 3
b <- 1000
N1 <- 1000
S1 <- 4
N2 <- 10
S2 <- 0
p <- seq(from=0,to=1,length=10000)
p1 <- dbeta(p,S1+a,N1-S1+b)
p2 <- dbeta(p,S2+a,N2-S2+b)
matplot(p[1:50],cbind(p1,p2)[1:50,],type="l")
(S1+a)/(N1+b)
(S2+a)/(N2+b)

• ついでRmd用テキストファイル

副作用のためのベータ分布
========================================================

>アレルギーはall or nothingで、何も起きないときにはなにも起こらないのですが、それはnext patientが起こさないことを意味しないので、別の安全な橋があるときはそっち渡ったほうがいいと思います。

という意見があったとします。

統計屋はこんな風に考えています(走り書きなので、細かいところは間違っているかもしれません)。

-Na回利用されて、Ka回の事件が起きた。そのうち0<=Sa<=Ka<=Na 回が報告された。
-Nb回利用されて、Kb回の事件が起きた。そのうち0<=Sb<=Kb<=Nb 回が報告された。
-Nc回利用されて、Kc=0回の事件が起きた。そのうち0=Sc<=Kc=0<=Nc 回が報告された。

というときに「nothingな安全な橋」とは、Ki=0が知りえないことが現実であるのでSi=0を持って代用すると思います。

今、(N,S)のデータの組みがいくつもあったとき、それに基づいてK/Nの値がいくつだろう、と思いを巡らします。

-(N1,S1)=(1000,4)
-(N2,S2)=(10,0)

の２つで考えたときに、

>K/N=0の確率はK1/N1=0の可能性はないのに対して、

>K2/N2=0の可能性はあるので、

"nothing"である可能性で選ぶなら、(N2,S2)=(10,0)の道を選ぶと思います。

そうは言っても(N2,S2)を見て、K2/N2=0.3のようなものすごく高い確率も否定はできないです。

>じゃあ、(N1,S1),(N2,S2)とを見たとき、どちらが、次の利用のときに事件が起きる確率は高いのか…

と考えたとき、ベータ分布を利用することがあります。
K/Nの値は0-1の値を取りますが、その値が、大小どんな値になるか全く見当もつかないときには:

>(N1,S1)の場合の「次の利用のときの事件発生率」

> 0.00499002

>(N2,S2)の場合の「次の利用のときの事件発生率」

> 0.08333333

そうはいっても副作用のようにめったにおきるものではない、と思う根拠(治験を通っている…とか)があるときには、K/Nの値として0-1の値のどれもが平等に起きるとは思えないです。かなり小さい値だろうと思います。そんなときは、「リスク想定」の下駄を履かせます(事前分布を与える)。

0.3%くらい？？とあたりをつければ
>(N1,S1)の場合の「次の利用のときの事件発生率」

> 0.0035

>(N2,S2)の場合の「次の利用のときの事件発生率」

> 0.002970297
 
と、逆転しました。

図は下駄を履かせた場合。黒は(N1,S1)、赤は(N2,S2)で、K/Nが低確率(0-0.5%)の範囲を描いています。
カーブの下面積を横軸の値で重み付き積分したのが「期待値」です。

あくまでも「ベータ分布」に準えたベイズ的な判断の話です。
医学判断はベイズ的なことが非常に多いので、その基礎である二項分布とベータ分布との関係は丁寧に勉強すると得すると思います。


```{r}
N1 <- 1000
S1 <- 4
N2 <- 10
S2 <- 0
# K/Nの値の想定値をいろいろとる
p <- seq(from=0,to=1,length=10000)
# ベータ分布に準えて、(N,S)の値をあてはめる
p1 <- dbeta(p,S1+1,N1-S1+1)
p2 <- dbeta(p,S2+1,N2-S2+1)
matplot(p[1:50],cbind(p1,p2)[1:50,],type="l")
# それぞれの次の利用時の事件の発生確率
(S1+1)/(N1+2)
(S2+1)/(N2+2)
```

```{r fig.width=7, fig.height=6}
# riskの下駄
# K/Nの値は0-1でどれかわからない、ということはないでしょう、かなり0に近いでしょう、というときは下駄を履かせます
a <- 3
b <- 1000
N1 <- 1000
S1 <- 4
N2 <- 10
S2 <- 0
p <- seq(from=0,to=1,length=10000)
p1 <- dbeta(p,S1+a,N1-S1+b)
p2 <- dbeta(p,S2+a,N2-S2+b)
matplot(p[1:50],cbind(p1,p2)[1:50,],type="l")
(S1+a)/(N1+b)
(S2+a)/(N2+b)
```