月: 2025年11月

#####Bayesian meta-analysis of Hazard Raio (HR)######

library(rstan)
library(bayesplot)
library(dplyr)
library(HDInterval)
library(forestplot)

# Get data via clipboard.
exdato = read.delim("clipboard",sep="\t",header=TRUE)
labem = exdato[,"label"]
em = labem[6]
labe_study = labem[1]
labe_cont = labem[2]
labe_int = labem[3]
labe_outc = labem[4]
labe_em = labem[5]

exdat=na.omit(exdato)
# Number of studies
K <- nrow(exdat)

# データの準備
if(colnames(exdat)[1]=="author"){
study = exdat$author
}
if(colnames(exdat)[1]=="study"){
study = exdat$study
}
nc=exdat$nc
nt=exdat$nt
meta_data <- list(
  K = K,
  yi = exdat$yi,
  sei = exdat$sei
)

#Stan model code
stan_model_code <- "
data {
  int<lower=0> K;              // 研究数
  vector[K] yi;                // 各研究の logHR
  vector<lower=0>[K] sei;      // 各研究の標準誤差
}
parameters {
  real mu;                     // 全体平均 logHR
  real<lower=0> tau;           // 異質性（標準偏差）
  vector[K] theta;             // 各研究の真の効果
}
model {
  theta ~ normal(mu, tau);     // 真の効果の分布
  yi ~ normal(theta, sei);     // 観測値の分布
}
generated quantities {
  real tau_sq = tau^2;         // 異質性の分散
  real I2 = tau_sq / (tau_sq + mean(sei .* sei)); // I²の推定
  real pred = normal_rng(mu, sqrt(tau_sq + mean(sei .* sei))); // 予測区間の一例
}"


# Stan モデルのコンパイルとサンプリング
fit <- stan(
  model_code = stan_model_code,
  data = meta_data,
  iter = 30000,
  warmup = 10000,
  chains = 4,
  seed = 123,
  control = list(adapt_delta = 0.98)
)

# 結果の表示
print(fit, pars = c("mu", "tau", "tau_sq", "I2", "pred"), probs = c(0.025, 0.5, 0.975))

# 事後分布の可視化
posterior <- extract(fit)
#mcmc_areas(as.data.frame(posterior), pars = c("mu", "tau", "I2", "pred"))

# Summary HR and 95% credible interval
mu_HR = exp(mean(posterior$mu))
mu_HR_lw = exp(quantile(posterior$mu,probs=0.025))
mu_HR_up = exp(quantile(posterior$mu,probs=0.975))
p_val_HR = round(2*min(mean(posterior$mu >0), mean(posterior$mu <0)), 6)	#Bayesian p-value for summary logRR
if(mean(posterior$mu)>0){
prob_direct_HR = round(mean(posterior$mu > 0),6)
label = "p(HR>1)="
}else{
prob_direct_HR = round(mean(posterior$mu < 0),6)
label = "p(HR<1)="
}
# Prediction Inerval and 95% credible interval
pi_HR_lw = exp(quantile(posterior$pred,probs=0.025))
pi_HR_up = exp(quantile(posterior$pred,probs=0.975))

# Tau-squared mode and 95% HDI (High Density Interval) of logHR
dens=density(posterior$tau_sq)
tau_squared=dens$x[which.max(dens$y)]	#mode
ci_tau_squared=hdi(posterior$tau_sq)	#95% HDI

# I-squared mode and 95% HDI (High Density Interval) of logHR
dens=density(posterior$I2)
I_squared=100*dens$x[which.max(dens$y)]		#mode
ci_I_squared=100*hdi(posterior$I2)		#95% HDI

# logHR of each study = theta[,i]
hr=rep(0,K)
hr_lw=rep(0,K)
hr_up=rep(0,K)
for(i in 1:K){
hr[i] = exp(mean(posterior$theta[,i]))
hr_lw[i] = exp(quantile(posterior$theta[,i],probs=0.025))
hr_up[i] = exp(quantile(posterior$theta[,i],probs=0.975))
}
# 各研究の効果推定値の重み%各研究の効果推定値の重み%:事後不確実性の程度（推定の安定性）
# 重みを格納するベクトルを初期化
weights <- numeric(K)
weight_percentages <- numeric(K)
for (i in 1:K) {
  # i番目の研究のlogRRのサンプリング値を取得
  theta_i_samples <- posterior$theta[, i]
  # そのサンプリングされた値の分散を計算
  # これを「実効的な逆分散」と考える
  variance_of_theta_i <- var(theta_i_samples)
  # 重みは分散の逆数
  weights[i] <- 1 / variance_of_theta_i
}
# 全重みの合計
total_weight <- sum(weights)
# 各研究の重みのパーセンテージ
weight_percentages <- (weights / total_weight) * 100
wpc = format(round(weight_percentages,digits=1), nsmall=1)
wp = weight_percentages/100
#Forest plot box sizes on weihts
wbox=c(NA,NA,(K/4)*sqrt(wp)/sum(sqrt(wp)),0.5,0,NA)	

#Forest plot by forestplot
#setting fs for cex
fs=1
if(K>20){fs=round((1-0.02*(K-20)),digits=1)}

m=c(NA,NA,hr,mu_HR,mu_HR,NA)
lw=c(NA,NA,hr_lw,mu_HR_lw,pi_HR_lw,NA)
up=c(NA,NA,hr_up,mu_HR_up,pi_HR_up,NA)

hete1=""
hete2=paste("I2=",round(I_squared, 2),"%",sep="")

hete3=paste("tau2=",format(round(tau_squared,digits=4),nsmall=4),sep="")

#hete2=paste("I2=",round(mean_I_squared_logRR, 2), "%\n"," (",round(ci_I_squared_logRR[1], 2)," ~ ",
#round(ci_I_squared_logRR[2], 2),")",sep="")

#hete3=paste("tau2=",format(round(mean_tau_squared_logRR,digits=4),nsmall=4), "\n",
#" (",format(round(ci_tau_squared_logRR[1],digits=4),nsmall=4)," ~ ",
#format(round(ci_tau_squared_logRR[2],digits=4),nsmall=4),")  ",sep="")

hete4=""
hete5=paste("p=",p_val_HR, sep="")
hete6=paste(label,prob_direct_HR,sep="")

au=study
sl=c(NA,toString(labe_study),as.vector(au),"Summary Estimate","Prediction Interval",NA)

spac=c("    ",NA,rep(NA,K),NA,NA,NA)

ncl=c(NA,"Number",nc,NA,NA,hete1)
rcl=c(NA,NA,rep(NA,K),NA,NA,hete2)

ntl=c(labe_int,"Number",nt,NA,NA,hete3)
rtl=spac

ml=c(labe_outc,labe_em,format(round(hr,digits=3),nsmall=3),format(round(mu_HR,digits=3),nsmall=3),NA,hete5)
ll=c(NA,"95% CI lower",format(round(hr_lw,digits=3),nsmall=3),format(round(mu_HR_lw,digits=3),nsmall=3),format(round(pi_HR_lw,digits=3),nsmall=3),hete6)
ul=c(NA,"95% CI upper",format(round(hr_up,digits=3),nsmall=3),format(round(mu_HR_up,digits=3),nsmall=3),format(round(pi_HR_up,digits=3),nsmall=3),NA)
wpcl=c(NA,"Weight(%)",wpc,100,NA,NA)
ll=as.vector(ll)
ul=as.vector(ul)
sum=c(TRUE,TRUE,rep(FALSE,K),TRUE,FALSE,FALSE)
zerov=1
labeltext=cbind(sl,ntl,rtl,ncl,rcl,spac,ml,ll,ul,wpcl)
hlines=list("3"=gpar(lwd=1,columns=1:11,col="grey"))
dev.new()
plot(forestplot(labeltext,mean=m,lower=lw,upper=up,is.summary=sum,graph.pos=7,
zero=zerov,hrzl_lines=hlines,xlab=toString(labe_em),txt_gp=fpTxtGp(ticks=gpar(cex=fs),
xlab=gpar(cex=fs),cex=fs),xticks.digits=2,vertices=TRUE,graphwidth=unit(50,"mm"),colgap=unit(3,"mm"),
boxsize=wbox,
lineheight="auto",xlog=TRUE,new_page=FALSE))
####
###Pint indiv estimates with CI and copy to clipboard.
nk=K+2
efs=rep("",nk)
efeci=rep("",nk)
for(i in 1:nk){
efs[i]=ml[i+2]
efeci[i]=paste(ll[i+2],"~",ul[i+2])
}

efestip=data.frame(cbind(efs,efeci))
efestip[nk,1]=""
print(efestip)
write.table(efestip,"clipboard",sep="\t",row.names=FALSE,col.names=FALSE)
print("The estimates of each study and the summary estimate are in the clipboard.")

#############

########################
library(metafor)

#Plot funnel plot.
if(em == "HR"){
yi = hr
vi=1/weights
sei = sqrt(1/weights)
mu = mu_HR
dev.new(width=7,height=7)
funnel(yi, sei = sei, level = c(95), refline = mu,xlab = "HR", ylab = "Standard Error")
}

#Asymmetry test with Egger and Begg's tests.
egger=regtest(yi, sei=sei,model="lm", ret.fit=FALSE)
begg=ranktest(yi, vi)

#Print the results to the console.
print("Egger's test:")
print(egger)
print("Begg's test")
print(begg)

#Add Begg and Egger to the plot.
fsfn=1
em=toString(exdat$label[6])
outyes=toString(exdat$label[4])
funmax=par("usr")[3]-par("usr")[4]
gyou=funmax/12
fxmax=par("usr")[2]-(par("usr")[2]-par("usr")[1])/40
text(fxmax,gyou*0.5,"Begg's test",pos=2,cex=fsfn)
kentau=toString(round(begg$tau,digits=3))
text(fxmax,gyou*1.2,paste("Kendall's tau=",kentau,sep=""),pos=2,cex=fsfn)
kenp=toString(round(begg$pval,digits=5))
text(fxmax,gyou*1.9,paste("p=",kenp,sep=""),pos=2,cex=fsfn)

text(fxmax,gyou*3.5,"Egger's test",pos=2,cex=fsfn)
tstat=toString(round(egger$zval,digits=3))
text(fxmax,gyou*4.2,paste("t statistic=",tstat,sep=""),pos=2,cex=fsfn)
tstap=toString(round(egger$pval,digits=5))
text(fxmax,gyou*5.1,paste("p=",tstap,sep=""),pos=2,cex=fsfn)
#####################################