Phân tích tổng hợp bằng BAYES

1 Giới thiệu và thí dụ minh họa

Phân tích tổng hợp (Meta-analysis) có ý nghĩa rất quan trọng đối với Y học thực chứng. Nhà nghiên cứu sẽ tổng hợp kết quả từ nhiều nghiên cứu riêng lẻ (thường là thử nghiệm lâm sàng) có cùng loại outcome , thí dụ giá trị trung bình, Odds-ratio, RR, incidence rate ratio… đã được công bố trong y văn, từ đó ước tính chính xác hơn hiệu ứng thực sự sau khi hiệu chỉnh cho sai biệt do ngẫu nhiên cũng như đánh giá mức độ đồng nhất giữa các nghiên cứu đơn lẻ.

Bài thực hành hôm nay chỉ giới hạn cho một dạng đặc biệt trong phân tích tổng hợp, đó là sử dụng mô hình hồi quy đa cấp (hierarchical model) hoặc mô hình có hiệu ứng ngẫu nhiên (random effect) để ước tính kích thước hiệu ứng trung bình. Nhưng đặc biệt hơn, chúng ta sẽ thực hiện mô hình này theo phương pháp Bayes, chứ không dùng phương pháp Maximum Likelihood estimation (ước tính hợp lý cực đại) của phái frequentist như truyền thống.

Thí dụ minh họa trong bài là một nghiên cứu của Graham Colditz đặng trên tờ JAMA năm 1994, nhằm ước lượng hiệu quả của vaccine BCG trong việc ngừa bệnh Lao. Trong nghiên cứu này, các tác giả phân tích tổng hợp 14 thử nghiệm lâm sàng từ năm 1949 đến 1980. Chỉ số thống kê chính được khảo sát là nguy cơ tương đối (Relative risk : RR) cho bệnh Lao dù có chủng ngừa, từ đó cho phép suy ra hiệu quả bảo vệ (Protective effect : PE= 1-RR).

Đây là ảnh chụp bài báo gốc, trong bảng ghi lại toàn bộ dữ liệu thô trích xuất từ 13 nghiên cứu, dựa vào đó một mô hình hiệu ứng ngẫu nhiên (random effect) đã được dựng để ước tính RR tổng hợp và khoảng tin cậy của chỉ số này.

Dữ liệu thô này còn có thể trích xuất từ package metafor trong R :

library(tidyverse)

library(metafor)

dat <- get(data(dat.bcg))%>%.[,1:7]

head(dat)%>%knitr::kable()

trial	author	year	tpos	tneg	cpos	cneg
1	Aronson	1948	4	119	11	128
2	Ferguson & Simes	1949	6	300	29	274
3	Rosenthal et al	1960	3	228	11	209
4	Hart & Sutherland	1977	62	13536	248	12619
5	Frimodt-Moller et al	1973	33	5036	47	5761
6	Stein & Aronson	1953	180	1361	372	1079

2 Chuẩn bị dữ liệu và mô tả

Như ta thấy, dữ liệu thô được trình bày dưới dạng bảng chéo 2x2 (4 cột cho tần suất Có/không mắc bệnh Lao ở 2 phân nhóm Có/không tiêm chủng BCG)

t tương ứng với Treatment (có chủng ngừa), c tương ứng với Control (không chủng ngừa), pos và neg tương ứng với kết quả có /không mắc bệnh lao.

Thí dụ tpos tương ứng với những trường hợp có chủng ngừa BCG nhưng vẫn mắc bệnh lao

Từ những bảng chéo này, chúng ta có thể ước tính Relative risk (RR) và Protective effect dễ dàng như sau:

slab = paste(dat$author,dat$year,sep="-")

d1=dat$tpos
n1=dat$tpos+dat$tneg
d0=dat$cpos
n0=dat$cpos+dat$cneg
p1=d1/n1
p0=d0/n0

rr=p1/p0
pe=1-rr

cbind(Studies=slab, Relative_risk=rr,Protective_effect=pe)%>%knitr::kable()

Studies	Relative_risk	Protective_effect
Aronson-1948	0.410938654841094	0.589061345158906
Ferguson & Simes-1949	0.204868154158215	0.795131845841785
Rosenthal et al-1960	0.25974025974026	0.74025974025974
Hart & Sutherland-1977	0.236560523606413	0.763439476393587
Frimodt-Moller et al-1973	0.804489533795327	0.195510466204673
Stein & Aronson-1953	0.455611144836826	0.544388855163174
Vandiviere et al-1973	0.197721021611002	0.802278978388998
TPT Madras-1980	1.01202404809619	-0.0120240480961924
Coetzee & Berjak-1968	0.625366345142315	0.374633654857685
Rosenthal et al-1961	0.253765465303927	0.746234534696073
Comstock et al-1974	0.712226836059195	0.287773163940805
Comstock & Webster-1969	1.56191619962637	-0.561916199626368
Comstock et al-1976	0.982835076874145	0.0171649231258549

Nhưng ta không thể dùng trực tiếp RR cho phân tích tổng hợp được, mà phải hoán chuyển RR thành Log(RR).

Việc hoán chuyển RR bằng hàm logarit nhắm đến 2 mục tiêu, thứ nhất, để biến kết quả (outcome) mới là LogRR có thể được xấp xỉ bằng phân bố Gaussian, điều này rất quan trọng khi áp dụng phương pháp Bayes sau này ; thứ hai, với giả định LogRR phân bố chuẩn, ta có thể tính SE hay SD của logRR theo công thức :

\[\sigma^2_j \approx \frac{1 - p_{1j}}{n_{1j}p_{1j}} + \frac{1 - p_{0j}}{n_{0j}p_{0j}}\] Đây là kết quả:

lrr=log(rr)

lse <- sqrt((1 - p1)/(p1 * n1)+(1-p0)/(p0*n0))

ve<-lse^2

lower <- exp(lrr-qnorm(.975)*lse)
upper <- exp(lrr+qnorm(.975)*lse)

cbind(Studies=slab,LogRR=lrr,LRR_SE=lse)%>%knitr::kable()

Studies	LogRR	LRR_SE
Aronson-1948	-0.889311333920205	0.570600354892951
Ferguson & Simes-1949	-1.58538865720143	0.441113501718258
Rosenthal et al-1960	-1.34807314829969	0.644490469571091
Hart & Sutherland-1977	-1.44155119002131	0.141456819921302
Frimodt-Moller et al-1973	-0.217547322211296	0.226296646395016
Stein & Aronson-1953	-0.786115585818864	0.0831000508778951
Vandiviere et al-1973	-1.62089822359839	0.472247019654243
TPT Madras-1980	0.0119523335238412	0.0629410779842364
Coetzee & Berjak-1968	-0.469417648738149	0.237558856840255
Rosenthal et al-1961	-1.37134480347278	0.270231000466321
Comstock et al-1974	-0.33935882833839	0.111410116109623
Comstock & Webster-1969	0.445913400571378	0.729729981020482
Comstock et al-1976	-0.0173139482168795	0.267216503390016

Package metafor là một công cụ rất tốt để làm meta-analysis trong R, ta sẽ dùng nó để tóm tắt dữ liệu thô bằng forest plot.

ta có thể trình bày Protective effect

library("metafor")

pPE <- forest(x = 1-rr, 
            ci.lb = 1-upper,
            ci.ub = 1-lower,
            slab = paste(dat$author,
                         dat$year,sep="-"),
                         refline=0)

text(min(pPE$xlim), 0.9 * max(pPE$ylim), "Study and Year", pos = 4, font = 2)
text(max(pPE$xlim), 0.9 * max(pPE$ylim), "Protective effect [97.5% CI]", pos = 2, font = 2)

Hoặc RR

pRR <- forest(x = rr, 
              ci.lb = lower,
              ci.ub = upper,
              slab = paste(dat$author,
                           dat$year,sep="-"),
              refline=1)

text(min(pRR$xlim), 0.9 * max(pRR$ylim), "Study and Year", pos = 4, font = 2)
text(max(pRR$xlim), 0.9 * max(pRR$ylim), "Relative Risk [97.5% CI]", pos = 2, font = 2)

3 Mô hình đa cấp

Phân tích tổng hợp có thể được mô hình hóa bằng một mô hình hỗn hợp/đa cấp (hierarchical model), trong đó biến số kết quả, thí dụ LogRR cho mỗi nguyên cứu gốc j được giả định có phân bố chuẩn, với 2 tham số : trung bình Thetaj và độ lệch chuẩn Sigma (hay phương sai : sigma^2)

\[LRR_j \sim N(\theta_j, \sigma_j)\]

Trong đó sigma được giả định có giá trị xác định (vì những thử nghiệm lâm sàng nguyên thủy được thiết kế với cỡ mẫu lớn), nhưng Trung bình Thetaj thì bất định và được xét như hiệu ứng ngẫu nhiên (random effect) trong mô hình đa cấp mà ta sẽ dựng

Giá trị random effect Thetaj này được giả định có phân bố chuẩn, với tham số vị trí trung tâm là µ và sai số ngẫu nhiên là Tau.

\[\theta_j \sim N(\mu, \tau)\]

Hay cũng có thể viết:

\[\theta_j = \mu + \tau * \eta_j\]

Như vậy, µ đại diện cho logRR trung bình, cũng là hằng số (intercept) trong mô hình và chính là mục tiêu mà ta cần xác định, còn tích số giữa Tau và hiệu ứng ngẫu nhiên eta_j sẽ quyết định giá trị Thetaj gần hay xa so với hằng số µ.

4 Giải pháp MLE

Ta có thể dựng mô hình này trong metafor, sử dụng biến kết quả là LRR và phương sai của nó (ve=lse^2)

mod <- rma(lrr, ve)

summary(mod)

## 
## Random-Effects Model (k = 13; tau^2 estimator: REML)
## 
##   logLik  deviance       AIC       BIC      AICc  
## -12.2024   24.4047   28.4047   29.3746   29.7381  
## 
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value):      0.5597
## I^2 (total heterogeneity / total variability):   92.22%
## H^2 (total variability / sampling variability):  12.86
## 
## Test for Heterogeneity: 
## Q(df = 12) = 152.2330, p-val < .0001
## 
## Model Results:
## 
## estimate      se     zval    pval    ci.lb    ci.ub     
##  -0.7145  0.1798  -3.9744  <.0001  -1.0669  -0.3622  ***
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Kích thước hiệu ứng trung bình có thể tính cho cả RR và PE:

# RR
cbind(exp(mod$b),LB=exp(mod$ci.lb), UB=exp(mod$ci.ub))%>%knitr::kable()

		LB	UB
intrcpt	0.4894209	0.3440743	0.6961661

# Pef

cbind(1-exp(mod$b), LB=1-exp(mod$ci.ub),UB=1-exp(mod$ci.lb)) %>%knitr::kable()

		LB	UB
intrcpt	0.5105791	0.3038339	0.6559257

Một khi ước lượng được logRR từ mô hình, ta có thể dễ dàng phục hồi RR bằng exp(LogRR ), hoặc hiệu năng bảo vệ của BCG: Pef = 1-exp(LogRR).

5 Dựng mô hình BAYES bằng STAN

Ta viết một mô hình đơn giản như mô tả bằng STAN cod, sau đó cung cấp dữ liệu và lấy mẫu MCMC cho mô hình, chú ý là mô hình sử dụng lse và lrr

library("rstan")
set.seed(606)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

J <- nrow(dat)
stan.dat <- list(J = J, y = lrr, sigma = lse)

stan_code="
data {
  int<lower=0> J; // number of trials 
  real y[J]; // estimated log relative risk
  real<lower=0> sigma[J]; // se of log relative risk
}

parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}
"

Mô hình được lấy mẫu với 4 chuỗi, mỗi chuỗi 1000 lượt

fit <- stan(model_code = stan_code,
            data = stan.dat, 
            iter = 1500, 
            warmup = 500,
            control=list(adapt_delta=0.95),
            refresh = 0,
            chains = 4)

## In file included from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file1d142b8c31aa.cpp:8:
## C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"

Đây là bảng tóm tắt kết quả mô hình: Ta có Mu chính là Intercept, Tau là một hằng số trong khi tính hỗn tạp ngẫu nhiên của dữ liệu được ước tính từ Tau và 13 tham số etaj với j=1,2…13; cuối cùng, 13 giá trịThetaj = hiệu ứng logRR cá thể, ước tính cho 13 nghiên cứu gốc

library(broom)

tidyMCMC(fit,conf.int = TRUE, 
         conf.method = "HPDinterval")

##         term     estimate  std.error    conf.low   conf.high
## 1         mu -0.720848526 0.20808853 -1.15668263 -0.34315850
## 2        tau  0.642230441 0.18955901  0.33530760  1.02776004
## 3     eta[1] -0.158550400 0.70000980 -1.63003034  1.12814204
## 4     eta[2] -0.914730326 0.65931317 -2.26229055  0.31477161
## 5     eta[3] -0.474093856 0.73855951 -1.94168383  0.95655899
## 6     eta[4] -1.142694126 0.48214557 -2.09383653 -0.21549761
## 7     eta[5]  0.716440681 0.46722093 -0.23571796  1.60594037
## 8     eta[6] -0.115083278 0.34639851 -0.79495687  0.55546540
## 9     eta[7] -0.897003770 0.65926184 -2.21814982  0.39584488
## 10    eta[8]  1.208712070 0.44096613  0.39524113  2.08326610
## 11    eta[9]  0.346669453 0.46203110 -0.52554533  1.25532342
## 12   eta[10] -0.885320375 0.51782321 -1.87469978  0.18989151
## 13   eta[11]  0.608617687 0.37785359 -0.11141955  1.33332797
## 14   eta[12]  0.775737601 0.78135611 -0.81274867  2.26519040
## 15   eta[13]  0.945005604 0.51076079 -0.03567848  1.94641065
## 16  theta[1] -0.816408629 0.42297702 -1.61583418  0.04159691
## 17  theta[2] -1.284941952 0.38453668 -2.04207165 -0.54665701
## 18  theta[3] -1.020760442 0.46538036 -1.88903960 -0.05086047
## 19  theta[4] -1.401589338 0.14324979 -1.67509740 -1.11658303
## 20  theta[5] -0.281840018 0.21406374 -0.68686006  0.13975201
## 21  theta[6] -0.784874693 0.08460213 -0.95149241 -0.62087004
## 22  theta[7] -1.281150552 0.39982063 -2.09594405 -0.51610489
## 23  theta[8]  0.003561596 0.06235807 -0.11818770  0.12319012
## 24  theta[9] -0.506156121 0.22155815 -0.92362680 -0.06177265
## 25 theta[10] -1.256211715 0.25396763 -1.75855646 -0.77915606
## 26 theta[11] -0.352061961 0.11042684 -0.56475723 -0.12916498
## 27 theta[12] -0.217398131 0.51199058 -1.20276928  0.83206019
## 28 theta[13] -0.141153900 0.25061542 -0.63856011  0.34046881

6 Dựng mô hình Bayes bằng brms

Một cách làm khác nếu bạn không muốn code thủ công, đó là sử dụng package brms, như sau:

Đây là một mô hình rất đơn giản (chỉ có hằng số, và hiệu ứng ngẫu nhiên dựa vào lse và study labels). Mô hình được compile vào STAN code và lấy mẫu 4000 lượt/4 chuỗi như trên.

Sau khi mô hình converge, ta có thể sử dụng hàm forest của package brmstools để vẽ forest plot.

library(brms)

bayesdat=data_frame(slab=paste(dat$author,
                               dat$year,sep="-"),
                    LRR=lrr,
                    LSE=lse)

brm_out <- brm(LRR | se(LSE)~ 1 + (1|slab),
               data = bayesdat, 
               iter = 1500, 
               warmup = 500, 
               chains = 4,
               control=list(adapt_delta=0.95),
               refresh = 0,
               cores = 4)

brmstools::forest(brm_out,
       show_data = TRUE,
       av_name = "log_RR",fill_ridge="gold")+
  geom_vline(xintercept=0,linetype=2)

So với mô hình thủ công, mô hình brms có nội dung đơn giản hơn, nó chỉ ước tính phân phối hậu nghiệm cho b_intercept = tham số Mu, sd_slab_intercept hay hằng số hỗn tạp Tau, và 13 tham số Thetaj cho 13 nghiên cứu gốc:

tidyMCMC(brm_out,conf.int = TRUE, 
         conf.method = "HPDinterval")

##                                           term    estimate std.error
## 1                                  b_Intercept -0.73280615 0.2001215
## 2                           sd_slab__Intercept  0.63265734 0.1890787
## 3               r_slab[Aronson-1948,Intercept] -0.07811040 0.4335132
## 4      r_slab[Coetzee.&.Berjak-1968,Intercept]  0.22349371 0.2856339
## 5    r_slab[Comstock.&.Webster-1969,Intercept]  0.50223381 0.5051678
## 6        r_slab[Comstock.et.al-1974,Intercept]  0.37968769 0.2259522
## 7        r_slab[Comstock.et.al-1976,Intercept]  0.59575666 0.3153126
## 8      r_slab[Ferguson.&.Simes-1949,Intercept] -0.54816463 0.3873810
## 9  r_slab[Frimodt-Moller.et.al-1973,Intercept]  0.44830124 0.2801872
## 10    r_slab[Hart.&.Sutherland-1977,Intercept] -0.66765066 0.2371173
## 11      r_slab[Rosenthal.et.al-1960,Intercept] -0.28447600 0.4701317
## 12      r_slab[Rosenthal.et.al-1961,Intercept] -0.51914539 0.3076752
## 13      r_slab[Stein.&.Aronson-1953,Intercept] -0.05271437 0.2128750
## 14           r_slab[TPT.Madras-1980,Intercept]  0.73617485 0.2090598
## 15     r_slab[Vandiviere.et.al-1973,Intercept] -0.53868301 0.4017779
##       conf.low   conf.high
## 1  -1.14156049 -0.35458563
## 2   0.31654947  0.99838646
## 3  -0.91677342  0.75402466
## 4  -0.31538001  0.78705379
## 5  -0.56740148  1.43598299
## 6  -0.05621595  0.82323072
## 7  -0.03275966  1.21319960
## 8  -1.32761255  0.19501340
## 9  -0.07156576  1.01404354
## 10 -1.16539369 -0.23937560
## 11 -1.21127566  0.65190077
## 12 -1.15749683  0.07167702
## 13 -0.48978478  0.34949767
## 14  0.33284635  1.15847996
## 15 -1.29645290  0.28204882

Ta có thể vẽ thủ công một biểu đồ Forest plot khác cho kết quả Protective effect hoặc RR thay vì logRR:

postdf=brm_out$fit%>%as.data.frame()%>%
  as_tibble()%>%
  mutate(.,Iteration=as.numeric(rep(c(1:1000),4)),
         Chain=as.factor(rep(c(1:4),each=1000)))

colnames(postdf)=c("Overall","SD",as.vector(as.vector(bayesdat$slab)),"Lp","Iter","Chain")

longpost=postdf%>%
  gather(Overall,`Aronson-1948`:`Comstock et al-1976`,
         key="Study",value="LRR")%>%
  mutate(PE=1-exp(.$LRR),RR=exp(.$LRR))

library(ggridges)

library(viridis)

  ggplot(data=longpost)+
  geom_density_ridges_gradient(aes(x=PE,y=reorder(Study,PE),fill=..x..),
                               alpha=0.5,
                               scale=1.5,
                               show.legend = F)+
  scale_fill_viridis(option="D")+
  scale_x_continuous(limits = c(-2,1.5),"Protective effect")+
  scale_y_discrete("Studies")+
  theme_bw()+
  geom_vline(xintercept = 0,linetype=2,size=1)+
  geom_vline(xintercept = 1-exp(median(postdf$Overall)),linetype=2,col="blue")

 ggplot(data=longpost)+
  geom_density_ridges_gradient(aes(x=RR,y=reorder(Study,RR),fill=..x..),
                               alpha=0.5,
                               scale=1.5,
                               show.legend = F)+
  scale_fill_viridis(option="C",direction = -1)+
  scale_x_continuous(limits = c(0,6),"RR")+
  scale_y_discrete("Studies")+
  theme_bw()+
  geom_vline(xintercept = 1,linetype=2,size=1)+
  geom_vline(xintercept = exp(median(postdf$Overall)),linetype=2,col="blue")

## Picking joint bandwidth of 0.0571

7 Kết luận

Phân tích tổng hợp chỉ là một dạng đặc biệt của mô hình tuyến tính tổng quát có chứa hiệu ứng ngẫu nhiên, trong đó hiệu ứng cần khảo sát nằm ở Intercept, trong khi hiệu ứng bất định do tính hỗn tạp từ dữ liệu gốc được ước lượng bằng một tham số hiệu ứng ngẫu nhiên.

Như ta thấy, phương pháp Bayes mang lại một ưu thế mà phân tích tổng hợp truyền thống không có, đó là kích thước hiệu ứng được mô tả bằng một phân bố hậu nghiệm liên tục thực sự, chứ không chỉ giới hạn ở vài điểm cục bộ.

Việc sử dụng mô hình còn cho phép thực hiện những phân tích tổng hợp phức tạp khác trong đó biến số kết quả có những phân bố khác, thí dụ Poisson hay Negative binomial cho incidence rate, ta cũng có thể đưa thêm vào mô hình nhiều hiệu ứng ngẫu nhiên khác ngoài biến định danh tên nghiên cứu slab, thí dụ điều kiện địa lý, quốc gia, chủng tộc, giới tính, phương pháp phân tích dữ liệu, thiết kế …