時間イベントアウトカムつまり生存分析のハザード比Hazard Ratio (HR)の統合値と95％確信区間Credible Intervalを、Stanを用いてベイジアンメタアナリシスで求め、Forest plotを作成するコードです。
データはExcelで図のように準備します。介入群の症例数、対照群の症例数、ハザード比HRの自然対数、その標準誤差の値が必要になります。labelは研究名のラベル、対照、介入のラベル、アウトカム、効果指標のタイプ、その略称です。
使用法は今まで紹介した、RR, OR, HR, MD, SMDのためのスクリプトと同じです。Rを起動しておき、エディターにコピーした下記コードを貼り付けて、最初の5つのライブラリの読み込みを実行しておきます。Excelに戻って、セルB3からG9までの範囲を選択して、コピー操作（Ctrl+C)を行って、Rに戻り、exdato = read.delim(“clipboard”,sep=”\t”,header=TRUE)のスクリプトを実行します。続いて、それ以下のスクリプトを全部実行させると、Forest plotが出力されます。効果指標の点推定値と95％確信区間、および予測区間を提示します。
Markov Chain Monte Carlo (MCMC)シミュレーションを行うので、少し時間がかかります。chains = 4, warmup = 10000, iter = 30000に設定しているので、計80000個がサンプリングされます。ウォームアップ、繰り返しの回数は必要に応じて変更してください。
Forest plotが出力された時点で、各研究の効果推定値、95％確信区間、統合値と95％確信区間、予測区間の値をクリップボードにコピーしていますので、Excel評価シートなどに貼り付けることができます。
ハザード比HRのメタアナリシス用コード
#####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.")
#############
Forest plotが出力された後、続けて以下のFunnel plot作成用のスクリプトを実行すると、Funnel plotが出力されます。
Funnel plot作成用のスクリプト
########################
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)
#####################################
Funnel plotおよび非対称性の検定であるBeggの検定、Eggerの検定を出版バイアスの評価に用いる場合、研究数が10以上（5以上という研究者もいる）必要とされています。