# 日本計量生物学会 2023年度計量生物セミナー
# 「階層ベイズモデルとその応用」
#  野間 久史　（統計数理研究所）
# 2023年12月14日

#  R example code for performing Bayesian analyses of the meta-analysis examples in Higgins and Spiegelhalter (Int J Epidelol. 31(1):96-104).

###

install.packages("R2OpenBUGS")			# Installation command
install.packages("coda")

###

library("R2OpenBUGS")

# あらかじめ、以下のページから、OpenBUGSをインストールしておく必要があります。
# OpenBUGSは、デフォルトで、バックグラウンド稼働となりますので、すべての操作をR上で完結することができます。
# https://www.mrc-bsu.cam.ac.uk/software/bugs/openbugs/


###

## データセット(1)： BMJのメタアナリシスにおける8試験のデータセット

rc <- c(2,23,7,1,8,9,3,118)
rm <- c(1,9,2,1,10,1,1,90)
nc <- c(36,135,200,46,148,56,23,1157)
nm <- c(40,135,200,48,150,59,25,1159)
N <- 8

madata8 <- list(rc=rc,rm=rm,nc=nc,nm=nm,N=N)

###

## データセット(2)： ISIS-4 を除く14試験のデータセット

rc <- c(2,23,7,1,8,9,3,118,1,11,7,13,8,17)
rm <- c(1,9,2,1,10,1,1,90,0,6,1,5,4,4)
nc <- c(36,135,200,46,148,56,23,1157,21,75,27,33,122,108)
nm <- c(40,135,200,48,150,59,25,1159,22,76,27,23,130,107)
N <- 14

madata14 <- list(rc=rc,rm=rm,nc=nc,nm=nm,N=N)

###

## データセット(3)： ISIS-4 を含む15試験のデータセット

rc <- c(2,23,7,1,8,9,3,118,1,11,7,13,8,17,2103)
rm <- c(1,9,2,1,10,1,1,90,0,6,1,5,4,4,2216)
nc <- c(36,135,200,46,148,56,23,1157,21,75,27,33,122,108,29039)
nm <- c(40,135,200,48,150,59,25,1159,22,76,27,23,130,107,29011)
N <- 15

madata15 <- list(rc=rc,rm=rm,nc=nc,nm=nm,N=N)

###

## BUGSのモデルコード（R上で、関数と定義して、R2OpenBUGSでそのまま使用することができます）

## (a) Random-effects model; reference priorによる解析

rma01 <- function(){

	for (i in 1:N) {
		rc[i] ~ dbin(pc[i],nc[i]);
		rm[i] ~ dbin(pm[i],nm[i]);
		phi[i] <- logit(pc[i]);
		logit(pm[i]) <- phi[i] + delta[i];
		delta[i] ~ dnorm(mu, precision);		# Random effects
		pc[i] ~ dunif(0,1);						# Prior for pc
	}

	delta.new ~ dnorm(mu, precision);			# Predicted effect
		
	mu ~ dnorm(0, 0.0001);						# Reference prior for mu
	tau ~ dunif(0, 100);

	precision <- 1/(tau*tau);
	tau.sq <- 1/precision;
	
	OR <- exp(mu);
	
	less0 <- min(delta.new, 0);					# Calculate probabilities
	prob0 <- 1 - equals(less0, 0);
	less10 <- min(delta.new, -0.1054);
	prob10 <- 1 - equals(less10, -0.1054);
		
}


## (b) Random-effects model; sceptical priorによる解析

rma02 <- function(){

	for (i in 1:N) {
		rc[i] ~ dbin(pc[i],nc[i]);
		rm[i] ~ dbin(pm[i],nm[i]);
		phi[i] <- logit(pc[i]);
		logit(pm[i]) <- phi[i] + delta[i];
		delta[i] ~ dnorm(mu, precision);		# Random effects
		pc[i] ~ dunif(0,1);						# Prior for pc
	}

	delta.new ~ dnorm(mu, precision);			# Predicted effect
		
	mu ~ dnorm(0, 32.69);						# Sceptical prior for mu
	tau ~ dunif(0, 100);

	precision <- 1/(tau*tau);
	tau.sq <- 1/precision;
	
	OR <- exp(mu);
	
	less0 <- min(delta.new, 0);					# Calculate probabilities
	prob0 <- 1 - equals(less0, 0);
	less10 <- min(delta.new, -0.1054);
	prob10 <- 1 - equals(less10, -0.1054);
		
}

		
###

## 初期値の設定

mainits <- function(){
	list(mu = 0, tau = 1)
}


## OpenBUGSによるMCMC: (1) Reference priorによる解析: 14 trials

rmaout.01 <- bugs(data=madata14,
	inits=mainits,
	parameters.to.save=c("OR","mu","tau","delta.new","prob0","prob10"), 
	model.file=rma01, 
	n.chains=1, 
	n.iter=22000, 
	n.burnin=2000)


## OpenBUGSによるMCMC: (2) Reference priorによる解析: 14 trials + ISIS-4

rmaout.02 <- bugs(data=madata15,
	inits=mainits,
	parameters.to.save=c("OR","mu","tau","delta.new","prob0","prob10"), 
	model.file=rma01, 
	n.chains=1, 
	n.iter=22000, 
	n.burnin=2000)


## OpenBUGSによるMCMC: (3) Sceptical priorによる解析: 14 trials

rmaout.03 <- bugs(data=madata14,
	inits=mainits,
	parameters.to.save=c("OR","mu","tau","delta.new","prob0","prob10"), 
	model.file=rma02, 
	n.chains=1, 
	n.iter=22000, 
	n.burnin=2000)


## OpenBUGSによるMCMC: (4) Sceptical priorによる解析: 14 trials + ISIS-4

rmaout.04 <- bugs(data=madata15,
	inits=mainits,
	parameters.to.save=c("OR","mu","tau","delta.new","prob0","prob10"), 
	model.file=rma02, 
	n.chains=1, 
	n.iter=22000, 
	n.burnin=2000)



## 結果のサマリーの出力
	
rmaout.01$summary
rmaout.02$summary
rmaout.03$summary
rmaout.04$summary



###

## coda packageのツールによる収束診断など

library("coda")

crmaout.03 <- bugs(data=madata14,
	inits=mainits,
	parameters.to.save=c("mu","tau","delta.new"), 
	model.file=rma02, 
	n.chains=5, 
	n.iter=5000, 
	n.burnin=2000,
	codaPkg=TRUE)				# "codaPkg=TRUE" の引数を入れておくと、codaとの連携ができるオブジェクトが出力されます。

ermaout.03 <- read.bugs(crmaout.03)

summary(ermaout.03)
HPDinterval(ermaout.03)			# 最高事後密度区間

batchSE(ermaout.03)				# バッチ標準誤差
autocorr.diag(ermaout.03)		# 自己相関係数による診断

plot(ermaout.03)				# Trace plot and density estimation

gelman.diag(ermaout.03)			# Gelman-Rubin診断
gelman.plot(ermaout.03)

geweke.diag(ermaout.03)			# Gewekeの収束診断
geweke.plot(ermaout.03)