Gamma-Poisson Shrinker (GPS) Program をRで実装してみよう~シグナル検出の手法~
この記事では、Gamma-Poisson Shrinker (GPS) Programについて説明しつつ、それをRで実装してみようと思います。
シグナル検出ってそもそもなんだろう
そもそも、シグナル検出とはなんでしょうか。シグナル検出は薬剤疫学という分野で使用されます。薬には副作用がつきもので、実際に市販された薬に関してさまざまな副作用(厳密には有害事象)の情報が報告されてきます。様々な副作用の情報の中には「ほかの薬剤と比べても薬剤Aは副作用が出すぎじゃないか・・・」と思われるような数字が出てきます。人間が見てすぐ判断できるような数字なら問題ありませんが(薬的には大問題ですが)、データ量が多すぎると人間にはうまく解釈できなかったりします。人間だけだと「一体どの副作用の情報から調査したらいいんだろう」と手をこまねてしまうことになります。そんな問題を解決するのが、「シグナル検出」というものです。
シグナル検出は「この薬剤のこの副作用が、ほかの薬剤より出すぎじゃないか。調査したほうがいいんじゃないかランキング」という順位付けを行うことを狙いとしています。
GPSの気持ちを考えよう
記号を置こう
それでは、統計的にシグナル検出の気持ちを考えます。シグナル検出で使用されるデータは分割表を用いて以下のように記述されます。有害事象の種類をN、薬剤の種類をMとします。データのセルの数はM*Nとなります。シグナル検出の役割はこの表が与えられたときに、「このって数字は大きすぎる!」と指摘するということになります。
名前 | 有害事象1 | 有害事象2 | ・・・ | 有害事象N | 合計 |
---|---|---|---|---|---|
薬剤1 | ・・・ | ||||
薬剤2 | ・・・ | ||||
・・・ | |||||
薬剤M | ・・・ | ||||
合計 | ・・・ |
次に、「どうなればが大きい」と判断できるでしょうか。それには比較対象が必要です。そこで、ベースラインを設定します。はデータのもとで、「たぶん、このぐらいの値だろうな」という期待値みたいなものです。
- がに比べてめちゃくちゃ大きい → はシグナルだ!
- がに比べてそんなに大きくない、小さい → はシグナルじゃない!
大きい、小さいの基準をはっきりさせるために、とします。こうすることで
- が●●より大きい! → はシグナルだ!
- が●●より小さい! → はシグナルじゃない!
と言い換えることができました。
統計モデルを考えよう。
さきほどまでで観測等々の記号を置くことができました。次に、統計モデルを置きます。初めに観測をポアソン分布、 の事前分布を混合ガンマ分布を置きます。
このの混合事前分布は、『データの中には「シグナルっぽい群」と「シグナルっぽくない群」の2種類がある』をイメージしています。これのおかげで柔軟なモデルを実現しています。なぜガンマ分布かというとポアソン分布の共役事前分布に対応しているからです。
事後分布を計算しよう。
今回共役事前分布を用いているので簡単な計算での事後分布を導出することができます。
ハイパーパラメータを計算しよう。
上記で導出された事後分布には、まだ未決定のパラメータが5つあります。これらは、初期値をと設定し最適化を行います。
シグナルを判定する。
の下側5%点を計算し、これをという基準とします。これがのときシグナルと判定するのが流儀のようです。
GPSを実装しました
set.seed(125) #install.packages("MCMCpack") library(MCMCpack) N<-10;#イベントの数 M<-5;#薬剤の種類 S<-2;#属性による層別。今回は男女の2種類としてみた。 #各薬剤のイベントの総数を(適当に)決める #男女の人数は異なる Men.drug.user.sum<-round(seq(50,10000,length=M)) Female.drug.user.sum<-round(seq(1000,5000,length=M)) #各薬剤ごとに各イベントが何件ずつ発生したかを作りたい #今回は男女でイベントの発生確率は同じとする dirichlet.parameter<-rep(1:5,length=N) event.prop.matrix <- rdirichlet(M, dirichlet.parameter ) Men.drug.event.matrix<-c() Female.drug.event.matrix<-c() for(kk in 1:M){ Men.drug.event.vec<-as.vector(rmultinom(1,Men.drug.user.sum[kk],event.prop.matrix[kk,])) Men.drug.event.matrix<-rbind(Men.drug.event.matrix,Men.drug.event.vec) Female.drug.event.vec<-as.vector(rmultinom(1,Female.drug.user.sum[kk],event.prop.matrix[kk,])) Female.drug.event.matrix<-rbind(Female.drug.event.matrix,Female.drug.event.vec) } #データのリストを作成する。 #今回はコードが動くかの検証なので、シグナルの確率は均一とする。 drug.event.list<-list(Men.drug.event.matrix,Female.drug.event.matrix) drug.event.list[[1]][2,1]<-drug.event.list[[1]][2,1]*5 drug.event.list[[1]][4,3]<-drug.event.list[[1]][4,3]*5 drug.event.list[[1]][2,6]<-drug.event.list[[1]][2,6]*5 drug.event.list[[1]][5,9]<-drug.event.list[[1]][5,9]*5 drug.event.list[[2]][2,1]<-drug.event.list[[2]][2,1]*5 drug.event.list[[2]][4,3]<-drug.event.list[[2]][4,3]*5 drug.event.list[[2]][2,6]<-drug.event.list[[2]][2,6]*5 drug.event.list[[2]][5,9]<-drug.event.list[[2]][5,9]*5 #層別の調整を行ったベースラインを計算する #n++s #n+js #ni+s n_ALL_ALL_s<-c(sum(drug.event.list[[1]]),sum(drug.event.list[[2]])) n_event_ALL_s<-list(apply(drug.event.list[[1]],2,sum),apply(drug.event.list[[2]],2,sum)) n_ALL_drug_s<-list(apply(drug.event.list[[1]],1,sum),apply(drug.event.list[[2]],1,sum)) baseline.matrix<-matrix(0,M,N) for(ii in 1:N){ for(jj in 1:M){ baseline.matrix[jj,ii]<-n_event_ALL_s[[1]][ii]*n_ALL_drug_s[[1]][jj]/n_ALL_ALL_s[[1]]+ n_event_ALL_s[[2]][ii]*n_ALL_drug_s[[2]][jj]/n_ALL_ALL_s[[2]] } } obs.matrix<-(drug.event.list[[1]]+drug.event.list[[2]]) obs.vec<-as.vector(obs.matrix) base.vec<-as.vector(baseline.matrix) Lf1<-function(nn,aa,bb,ee){ ans<-(1+bb/ee)^(-nn)*(1+ee/bb)^(-aa)*exp(lgamma(aa+nn)-lgamma(aa)-lgamma(nn+1)) return(ans) } Lf2<-function(nn,aa1,aa2,bb1,bb2,ee,pp){ ans<-sum(log(pp*Lf1(nn,aa1,bb1,ee)+(1-pp)*Lf1(nn,aa2,bb2,ee))) return(ans) } #ハイパーパラメータを決定するために最適化を行う aa1.opt<-aa2.opt<-c() bb1.opt<-bb2.opt<-c() pp.opt<-c() aa1.opt[1]<-0.2; bb1.opt[1]<-0.1 aa2.opt[1]<-2 ;bb2.opt[1]<-4 pp.opt[1]<-1/3 #aa1の最適化 tt<-1 for(tt in 1:100){ opt <- optimize(Lf2, interval = c(0, 50),maximum = TRUE, nn = obs.vec, #aa1=aa1.opt[tt], bb1=bb1.opt[tt], aa2=aa2.opt[tt], bb2=bb1.opt[tt], ee=base.vec, pp=pp.opt[tt]) aa1.opt<-c(aa1.opt,opt$maximum) #bb1の最適化 opt <- optimize(Lf2, interval = c(0, 50),maximum = TRUE, nn = obs.vec, aa1=aa1.opt[tt+1], #bb1=bb1.opt[tt], aa2=aa2.opt[tt], bb2=bb1.opt[tt], ee=base.vec, pp=pp.opt[tt]) bb1.opt<-c(bb1.opt,opt$maximum) #aa2の最適化 opt <- optimize(Lf2, interval = c(0, 50),maximum = TRUE, nn = obs.vec, aa1=aa1.opt[tt+1], bb1=bb1.opt[tt+1], #aa2=aa2.opt[tt], bb2=bb1.opt[tt], ee=base.vec, pp=pp.opt[tt]) aa2.opt<-c(aa2.opt,opt$maximum) #bb2の最適化 opt <- optimize(Lf2, interval = c(0, 50),maximum = TRUE, nn = obs.vec, aa1=aa1.opt[tt+1], bb1=bb1.opt[tt+1], aa2=aa2.opt[tt+1], #bb2=bb1.opt[tt], ee=base.vec, pp=pp.opt[tt]) bb2.opt<-c(bb2.opt,opt$maximum) #ppの最適化 opt <- optimize(Lf2, interval = c(0, 1),maximum = TRUE, nn = obs.vec, aa1=aa1.opt[tt+1], bb1=bb1.opt[tt+1], aa2=aa2.opt[tt+1], bb2=bb1.opt[tt+1], ee=base.vec#, #pp=pp.opt[tt] ) pp.opt<-c(pp.opt,opt$maximum) } len<-length(aa1.opt) optimized.aa1<-aa1.opt[len] optimized.bb1<-bb1.opt[len] optimized.aa2<-aa2.opt[len] optimized.bb2<-bb2.opt[len] optimized.pp<-pp.opt[len] post.f1<-Lf1(obs.vec,optimized.aa1,optimized.bb1,base.vec) post.f2<-Lf1(obs.vec,optimized.aa2,optimized.bb2,base.vec) Qn<-post.f1/(post.f1+post.f2) Elambda<-Qn*(optimized.aa1+obs.vec)/(optimized.bb1+base.vec)+(1-Qn)*(optimized.aa2+obs.vec)/(optimized.bb2+base.vec) EB05.vec<-c() EB05.x<-seq(0,5,by=0.001) for(kk in 1:length(obs.vec)){ EB05.est<-Qn[kk]*pgamma(EB05.x,optimized.aa1+obs.vec[kk],rate=optimized.bb1+base.vec[kk])+ (1-Qn[kk])*pgamma(EB05.x,optimized.aa2+obs.vec[kk],rate=optimized.bb2+base.vec[kk]) EB05.vec[kk]<-EB05.x[which((abs(EB05.est-0.05)==min(abs(EB05.est-0.05))))] } matrix(EB05.vec,M,N) signal.vec<-(c(2,4,2,6))+5*(c(1,3,6,9)-1) plot((1:50)[-signal.vec],EB05.vec[-signal.vec],ylim=c(0,5),xlim=c(0,50),col=1,xlab="",ylab="");par(new=T) plot((1:50)[signal.vec],EB05.vec[signal.vec],ylim=c(0,5),xlim=c(0,50),col=2,xlab="",ylab="");
結果を考えよう
今回のテストデータでは(i,j)=(2,1)(4,3)(2,6)(5,9)のデータを5倍することでシグナルと疑似的にしました(赤色の点です)。EP05の結果は以下のようになりました。
赤色が、シグナルとなるデータです。4個のうち3個はシグナルと判定されていますが、残る1個だけが2未満、すなわちシグナルではないと判定されています。
この結果は、シグナル検出という手法を用いてもすべてのシグナルを検知できるわけではないことを示しています。
今日のところは、なぜ一か所だけうまくいかなかったかの理由はわかりませんでした。
(実際、どの程度の確率で検出できない例が存在するかは、テストデータの組み方によっても異なるので先行研究の論文を見てもよくはわからない)
また、EB05はをベイズ解釈し補正したものと考えられます。そこでをplotしてみると EP05は直接を計算するよりも保守的な結果になっていることが確認できます。明らかに大きなもの以外はシグナルとみなさないようにする傾向があるのかもしれません。もう少し数値的シミュレーションで性質の検証が必要かもしれません。