Gamma-Poisson Shrinker (GPS) Program をRで実装してみよう~シグナル検出の手法~

この記事では、Gamma-Poisson Shrinker (GPS) Programについて説明しつつ、それをRで実装してみようと思います。

GPSの肝

  • 相対リスク(Relative Risk)にベイズ的解釈を適用し、シグナル郡と非シグナル郡の混合分布を作っている
  • 観測に層別を用いている(BCPNNやその他の手法は層別してないっぽい)

シグナル検出ってそもそもなんだろう

そもそも、シグナル検出とはなんでしょうか。シグナル検出は薬剤疫学という分野で使用されます。薬には副作用がつきもので、実際に市販された薬に関してさまざまな副作用(厳密には有害事象)の情報が報告されてきます。様々な副作用の情報の中には「ほかの薬剤と比べても薬剤Aは副作用が出すぎじゃないか・・・」と思われるような数字が出てきます。人間が見てすぐ判断できるような数字なら問題ありませんが(薬的には大問題ですが)、データ量が多すぎると人間にはうまく解釈できなかったりします。人間だけだと「一体どの副作用の情報から調査したらいいんだろう」と手をこまねてしまうことになります。そんな問題を解決するのが、「シグナル検出」というものです。

シグナル検出は「この薬剤のこの副作用が、ほかの薬剤より出すぎじゃないか。調査したほうがいいんじゃないかランキング」という順位付けを行うことを狙いとしています。

GPSの気持ちを考えよう

記号を置こう

それでは、統計的にシグナル検出の気持ちを考えます。シグナル検出で使用されるデータは分割表を用いて以下のように記述されます。有害事象の種類をN、薬剤の種類をMとします。データのセルの数はM*Nとなります。シグナル検出の役割はこの表が与えられたときに、「このn_{ij}って数字は大きすぎる!」と指摘するということになります。

名前 有害事象1 有害事象2 ・・・ 有害事象N 合計
薬剤1 n_{11} n_{12} ・・・ n_{1N} n_{1+}
薬剤2 n_{21} n_{22} ・・・ n_{2N} n_{2+}
\vdots \vdots \vdots ・・・ \vdots \vdots
薬剤M n_{M1} n_{M2} ・・・ n_{MN} n_{M+}
合計 n_{+1} n_{+2} ・・・ n_{+N} n_{++}

次に、「どうなればn_{ij}が大きい」と判断できるでしょうか。それには比較対象が必要です。そこで、ベースラインE_{ij}を設定します。E_{ij}はデータのもとで、「たぶん、このぐらいの値だろうな」という期待値みたいなものです。

  • n_{ij}E_{ij}に比べてめちゃくちゃ大きい → n_{ij}はシグナルだ!
  • n_{ij}E_{ij}に比べてそんなに大きくない、小さい → n_{ij}はシグナルじゃない!

大きい、小さいの基準をはっきりさせるために、\lambda_{ij}=n_{ij}/E_{ij}とします。こうすることで

  • \lambda_{ij}が●●より大きい! → n_{ij}はシグナルだ!
  • \lambda_{ij}が●●より小さい! → n_{ij}はシグナルじゃない!

と言い換えることができました。

ベースライン

ベースラインは層別を用いてGPSでは以下のように定義されます。なぜこのように記述するかは別記事で言及いたします。
E_{ij}=\sum_s \frac{n_{i+s}n_{+js}}{n_{++s}}


統計モデルを考えよう。

さきほどまでで観測等々の記号を置くことができました。次に、統計モデルを置きます。初めに観測をポアソン分布、 \lambda_{ij}の事前分布を混合ガンマ分布を置きます。


 f(n_{ij}| \lambda_{ij},E_{ij}) \sim {\rm Poi}(n_{ij}|\lambda E_{ij}),\quad (i=1,\ldots,N, j=1,\ldots,M)
\pi(\lambda_{ij}|\alpha_1,\beta_1,\alpha_2,\beta_2,p) \sim p{\rm Ga}(\lambda_{ij}|\alpha_1,\beta_1)+(1-p){\rm Ga}(\lambda_{ij}|\alpha_2,\beta_2)


この \lambda_{ij}の混合事前分布は、『データの中には「シグナルっぽい群」と「シグナルっぽくない群」の2種類がある』をイメージしています。これのおかげで柔軟なモデルを実現しています。なぜガンマ分布かというとポアソン分布の共役事前分布に対応しているからです。

事後分布を計算しよう。

今回共役事前分布を用いているので簡単な計算で \lambda_{ij}の事後分布を導出することができます。


 Q_n=\frac{P*NBi(n_{ij}|\alpha_1,\beta_1,E_{ij})}{P*NBi(n_{ij}|\alpha_1,\beta_1,E_{ij})+(1-P)*NBi(n_{ij}|\alpha_2,\beta_2,E_{ij})}
\pi( \lambda_{ij}|n)=Q_n {\rm Ga}(\lambda_{ij}|\alpha_1+n_{ij},\beta_1+E_{ij})+(1-p){\rm Ga}(\lambda_{ij}|\alpha_2+n_{ij},\beta_2+E_{ij})

ハイパーパラメータを計算しよう。

上記で導出された事後分布には、まだ未決定のパラメータ\alpha_1,\beta_1,\alpha_2,\beta_2,pが5つあります。これらは、初期値を\alpha_1=0.2,\beta_1=0.1,\alpha_2=2,\beta_2=4,p=1/3と設定し最適化を行います。

シグナルを判定する。

\pi( \lambda_{ij}|n)の下側5%点を計算し、これをEB05という基準とします。これがEB05 \geq 2のときシグナルと判定するのが流儀のようです。

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の結果は以下のようになりました。

f:id:AsaKiriSun:20190318213142p:plain
GPSによるシグナル検出の結果

赤色が、シグナルとなるデータです。4個のうち3個はシグナルと判定されていますが、残る1個だけが2未満、すなわちシグナルではないと判定されています。
この結果は、シグナル検出という手法を用いてもすべてのシグナルを検知できるわけではないことを示しています。
今日のところは、なぜ一か所だけうまくいかなかったかの理由はわかりませんでした。
(実際、どの程度の確率で検出できない例が存在するかは、テストデータの組み方によっても異なるので先行研究の論文を見てもよくはわからない)


また、EB05はn_{ij}/E_{ij}ベイズ解釈し補正したものと考えられます。そこでEB05/(n_{ij}/E_{ij})をplotしてみると

f:id:AsaKiriSun:20190318214906p:plain
EP05の縮約量
EP05は直接n_{ij}/E_{ij}を計算するよりも保守的な結果になっていることが確認できます。明らかに大きなもの以外はシグナルとみなさないようにする傾向があるのかもしれません。もう少し数値的シミュレーションで性質の検証が必要かもしれません。