Mô hình phân cấp phi tuyến tính BAYES

1 Đặt vấn đề

Có một số chuyên đề nghiên cứu đặc biệt trong ngành Y-Dược nhắm đến mục tiêu là khảo sát những đặc tính “động học”, thí dụ khả năng hấp thu/thải trừ một dược chất sau khi đưa vào cơ thể (Dược động học),sự chuyển biến nồng độ của một biomarker nội sinh trong luồng khí thở (volatile organic compound, VOC), hay sự phân bố của một chất khí trong phổi sau khi hít vào (khí động học), sự suy giảm một chức năng sinh lý do quá trình lão hóa …

Những bài toán này đều có đặc điểm chung là biến số kết quả được mô tả bởi một quy luật (hệ thống) động học phi tuyến tính phụ thuộc thời gian. Phương pháp hồi quy tuyến tính truyền thống không thể giải quyết được một hệ thống tương tự. Do đó, nghiên cứu sinh bắt buộc phải dùng những mô hình phi tuyến tính (non-linear models), thậm chí phức tạp hơn là mô hình phân cấp (hierarchical model) hoặc nhiều ngăn (compartimental model).

Trong bài thực hành hôm nay (thuộc Project Bayes for Vietnam), Nhi sẽ hướng dẫn các bạn dựng các mô hình phân cấp-phi tuyến tính theo trường phái Bayes, nhờ package “brms” của Paul Bürkner. Hai thí dụ minh họa sẽ được trình bày, một rất đơn giản dùng dữ liệu mô phỏng, một phức tạp hơn và có thực, liên quan đến mô hình Dược động học 1 ngăn (PK model).

2 Khởi động với một thí dụ đơn giản

Giả sử ta muốn ước lượng giá trị của một đại lượng y cho cá thể (i) tại một thời điểm ti bất kì trong khoảng thời gian T, ta cần một mô hình có dạng:

\[y_i \sim f(t_j ,\phi) + e_{ij} \quad , \quad 0\leq j \leq T\]

Trong đó Phi là một vector chứa các tham số hồi quy (parameters) của mô hình, còn e là sai số, j là số điểm khảo sát trong thí nghiệm, i là 1 cá thể được khảo sát (bệnh nhân, con vật). y là kết quả, được mô tả như biến số ngẫu nhiên và tuân theo quy luật phân bố giả định (thí dụ Gaussian). Mô hình được xác định bởi bản chất của hàm f với thời điểm tj là biến số (predictor).

Đặc điểm chung của các mô hình phi tuyến tính, đó là hàm f có thể chứa trong nó một hay nhiều hàm phi tuyến tính, thí dụ Exponential. Công thức mô hình với Phi và tj bị chi phối bởi hàm Exp này.

Trong thí dụ đơn giản nhất, giả định hàm f có nội dung như sau:

\[f = b_0 e^{-b1*t}\]

library(tidyverse)

my_theme <- function(base_size =8, base_family = "sans"){
  theme_bw(base_size = base_size, base_family = base_family) +
    theme(
      panel.grid.major = element_line(color = "gray"),
      panel.grid.minor = element_blank(),
      panel.background = element_rect(fill = NA),
      strip.background = element_rect(fill = "#3c0477", color = "#3c0477", size =0.5),
      strip.text = element_text(face = "bold", size = 8, color = "white"),
      legend.position = "bottom",
      legend.justification = "center",
      legend.background = element_blank()
    )
}

Ta sẽ mô phỏng một dữ liệu dựa theo hàm f đơn giản như trên, với giả định tham số b0=150 và b1=0.08 và thời gian T dài 60 giây

t<-c(1:60)
se=rnorm(60,2,0.5)
out <- rnorm(60,mean=150*exp(-0.08*t),se)

set.seed(1206)
df=data_frame(Time=as.numeric(t),Outcome=out)

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

Time	Outcome
1	142.67188
2	128.95137
3	114.52035
4	108.35965
5	100.55330
6	92.32339

ggplot(df,aes(Time,Outcome))+
  geom_smooth(col="red4",fill="red",alpha=0.5)+
  geom_point()+my_theme(10)

Hàm f này mô tả một sự suy giảm theo thời gian của Outcome, theo quy luật hàm exp(-b1*t).

Nhi sẽ dùng package brms để dựng một mô hình nhằm xác định phân bố hậu nghiệm của b0 và b1. Nhưng trước hết, Nhi phải cân nhắc về giả thuyết tiền định (prior) cho Outcome, b0 và b1.

Ta sẽ phác thảo mô hình trên giấy:

Đầu tiên ta giả định Outcome là 1 biến ngẫu nhiên có phân bố Gaussian, :

\[yi \sim N(\mu,\sigma)\] Sau đó ta phác thảo nội dung hàm f cho phép xác định tham số µ:

\[\mu_{j} \sim b0*exp(-b1*Time_j)\] Dựa vào biểu đồ, ta có thể phán đoán là b0 > 100 nhưng nhỏ hơn 200, do đó prior cho nó có thể là 1 phân bố Gaussian có trung bình = 150 và sd =50

Tương tự, ta dự đoán b1 có giá trị rất thấp, gần 0 nhưng không rõ là bao nhiêu, một prior gamma có thể thích hợp cho dự đoán này

curve(dnorm(x, 150,20), from = 0,300,main="b0 ~ N(150,20)",
      xlab="", ylab="", bty="n",col="red")

curve(dgamma(x,2,3), from = 0,5, main="b1~Gamma(2,3)",col="blue",
      xlab="", ylab="", bty="n")

library(brms)
library(rstan)

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

prior1 <- c(prior(normal(150,20), nlpar="beta0"),
  prior(gamma(2,3), nlpar = "beta1",lb=0))

prior1

##             prior class coef group resp dpar nlpar     bound
## 1 normal(150, 20)     b                      beta0          
## 2     gamma(2, 3)     b                      beta1 <lower=0>

Sau khi thiết lập giả định về prior, ta viết công thức cho mô hình phi tuyến tính trong brms. Lưu ý là kể từ 2017, brms cho phép viết công thức (formula) theo nhiều phong cách khác nhau, trước và trong khi code cho hàm brm. Nhi thích chuẩn bị công thức Trước khi viết tùy chỉnh cho hàm brm, vì cách làm việc theo công đoạn này cho phép hạn chế sai sót và tổ chức tốt hơn quy trình.

Có một điều thú vị, đó là brms cho phép khai báo thêm những tham số / hằng số khác ngoài những biến có sẵn trong dữ liệu trong khi viết công thức, điều này dẫn đến một khả năng thích nghi hết sức uyển chuyển của mô hình Bayes, cho phép ta mở rộng mô hình với nhiều thành phần phong phú, để mô tả một giả thuyết (hàm) bất kì, chứ không chỉ giới hạn ở các biến từ dữ liệu.

Thí dụ ta có thể viết công thức :

“Outcome ~ beta0exp(-beta1Time)”

Để khai báo 2 tham số b0, b1. Sau đó ta đưa công thức này vào hàm bf,đồng thời khai báo thêm nội dung rằng b0 và b1 sẽ được khảo sát như hằng số (Intercept) trong mô hình, sau cùng lưu ý đổi tùy chỉnh nl=TRUE để khẳng định đây là mô hình phi tuyến tính.

Sau khi đóng công thức, ta khai báo dữ liệu là df, và prior là prior_mod.

Cuối cùng là chế độ lấy mẫu cho chuỗi MCMC, Nhi cho sampler chạy 4 chuỗi ngắn 4x1000 lượt sau 500 lượt khởi động.

Ta có thể cho thi hành mô hình fit1

fml1<- Outcome ~ beta0*(exp(-beta1*Time))

fit1 <- brm(bf(fml1,
               beta0+beta1~1,
               nl = TRUE),
            data = df,
            family = brmsfamily("gaussian", link_sigma = "log"),
            prior = prior1,
            control = list(adapt_delta = 0.90, max_treedepth=20),
          seed = 1234, iter = 1500,chains = 4)

Mô hình converge rất nhanh vì đơn giản, ta có thể tóm tắt phân bố hậu nghiệm của 2 hằng số beta0 và beta1 cũng như sigma của Outcome

summary(fit1)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: Outcome ~ beta0 * (exp(-beta1 * Time)) 
##          beta0 ~ 1
##          beta1 ~ 1
##    Data: df (Number of observations: 60) 
## Samples: 4 chains, each with iter = 1500; warmup = 750; thin = 1;
##          total post-warmup samples = 3000
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## beta0_Intercept   151.84      1.43   149.05   154.59       1201 1.00
## beta1_Intercept     0.08      0.00     0.08     0.08       1415 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     2.33      0.22     1.95     2.79       1845 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fit1, ask = FALSE,theme=my_theme(8))

Như ta thấy, mô hình đã xác định beta0 =150.98 (từ 148.41 đến 153.74), beta1=0.08 và sigma = 2.21 (1.85 đến 2.67), kết quả này rất phù hợp với các tham số ta đã dùng để mô phỏng dữ liệu.

theme_set(my_theme())

marginal_effects(fit1)

3 Một thí dụ phức tạp hơn: Mô hình PK 1 ngăn

Tiếp theo chúng ta sẽ thử sức với một trường hợp có độ khó cao hơn, đó là một nghiên cứu dược động học của Theophylline : Vào năm 1994, Boeckmann, Sheiner và Beal đã công bố dữ liệu này từ kết quả thực nghiệm của BS. Robert Upton. Thuốc Theophylline được dùng trên 12 người theo đường uống, sau đó nồng độ Theophilline trong máu được khảo sát liên tục qua 11 thời điểm trong 25 giờ sau đó. Dữ liệu này sau đó được Pinheiro, J. C.và Bates, D. M. sử dụng để minh họa cho mô hình phi tuyến tính có hiệu ứng ngẫu nhiên vào năm 2000 (Mixed-effects Models in S and S-PLUS, Springer).

Ta tải dữ liệu và bắt đầu khảo sát nó

dat<-read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/datasets/Theoph.csv")


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

X	Subject	Wt	Dose	Time	conc
1	1	79.6	4.02	0.00	0.74
2	1	79.6	4.02	0.25	2.84
3	1	79.6	4.02	0.57	6.57
4	1	79.6	4.02	1.12	10.50
5	1	79.6	4.02	2.02	9.66
6	1	79.6	4.02	3.82	8.58

Dữ liệu gồm 5 biến, Subject để định danh cá thể, Wt chỉ cân năng (kg),Dose là liều thuốc uống (mg/kg), Time chỉ thời điểm khảo sát và conc là biến số kết quả cần nghiên cứu có đơn vị mg/L.

Đầu tiên, chúng ta sẽ tính lại Liều thuốc thực sự được đưa vào cơ thể bằng cách nhân cho thể trọng:

Khảo sát trực quan cho thấy noèéng độ Theophylline thay đổi theo quy luật như sau qua 11 thời điểm:

library(RColorBrewer)

coledges=colorRampPalette(c("red","blue","purple"))
mycol=(coledges(12))

sum_df=dat%>%group_by(Time)%>%
  summarise_at("conc",median)

dat%>%ggplot()+
  geom_line(aes(x=Time,y=conc,col=factor(Subject)),alpha=0.5)+
  geom_smooth(data=sum_df,aes(x=sum_df$Time,y=sum_df$conc),col="red4",fill="red",size=1,alpha=0.2)+scale_color_manual(values=mycol)

Để đơn giản hóa, chúng ta sẽ sử dụng mô hình dược động học 1 ngăn (One Compartment PK model), theo đó Lượng thuốc ban đầu (Dose) được giả định phân bố ngay tức khắc và đồng đều trong cơ thể - được xem như 1 ngăn chứa duy nhất với thể tích V (L), Nồng độ thuốc trong máu ở bất cứ thời điểm nào sau đó sẽ được xác định bởi 2 hằng số hấp thu (ka) và thải trừ (ke).

Giả định conc = biến ngẫu nhiên Y được xác định bằng phân bố Gaussian với 2 tham số Mu và sigma.

\[Y_{ij} \sim N(\mu_{ij},\sigma_{ij})\]

Cho mỗi cá thể i và thời điểm tj, giá trị Muij được ước lượng bằng 1 hàm f có nội dung :

\[\mu_{ij} \sim f(t_{j},\phi) = \frac{Dose*k_{a}}{V*(k_{a}-k_{e})}\left ( e^{-k_{e}*t_{j}} - e^{-k_{a}*t_{j}} \right )\]

Như vậy mục tiêu của chúng ta là xây dựng một mô hình cho phép xác định đồng thời giá trị của 2 hằng số ka và ke.

Nhi sẽ làm việc này theo từng bước:

1. Đầu tiên, Nhi đưa hàm f vào công thức mô hình phi tuyến tính. Hàm f này có nội dung gần giống như trên, tuy nhiên với giả định là Nhi chỉ muốn tính ka và ke (2 hằng số hấp thu/thải trừ) nhưng không quan tâm đến Volume, Nhi không nhân liều thuốc với cân nặng nhưng giữ nguyên giá trị Dose và nhân cho 1 biến tên là InV = 1/V.

fml <- conc~Dose*ka*InV/(ka-ke)*(exp(-ke*Time)-exp(-ka*Time))

fml

## conc ~ Dose * ka * InV/(ka - ke) * (exp(-ke * Time) - exp(-ka * 
##     Time))

2. Sau đó, Nhi đưa công thức này vào mô hình qua hàm bf, sau đó phát triển thêm 3 công thức nữa để ước tính InV ~ 1|Subject, ka và ke đều là hằng số (Intercept).

3. Sau đó Nhi áp dụng prior cho conc là N(Mu,sigma) với link function = log cho sigma và Identity cho Mu.

4. Tiếp theo, là bước quan trọng nhất, Nhi khai báo 3 prior cho 3 tham số InV, ka và ke. Có thể dự đoán là giá trị của chúng đều rất thấp, có thể nằm giữa 0 và 2; nên Nhi chọn prior là Gamma()

curve(dgamma(x, 1, 5), from = 0, 5, main="Gamma priors",
      xlab="", ylab="", bty="n",col="red")
curve(dgamma(x, 2,5), add = TRUE,col="blue")
curve(dgamma(x, 3,5), add=TRUE,col="green")
text(x = c(0.1,0.2,0.6), y = c(4,2,1.2), labels = c("InV", "ke", "ka"),col=c("red","blue","green"))

5. Sau cùng, Ta mô tả chế độ lấy mẫu cho sampler, bao gồm 2 chuỗi, mỗi chuỗi 1000 lượt

dat$Subject=as.factor(dat$Subject)

bpkmod <- brm(bf(fml,
             InV ~ 1|Subject,
             ka + ke ~ 1,
             nl = TRUE),
          data = dat, 
          family = brmsfamily("gaussian", link_sigma = "log"),
          prior = c(prior(gamma(1,5), nlpar = "InV", lb=0),
                    prior(gamma(3,5), nlpar = "ka", lb=0),
                    prior(gamma(2,5), nlpar = "ke", lb=0)),
          control = list(adapt_delta = 0.999, max_treedepth=30),
          seed = 1206, iter = 2000,warmup = 1000,chains = 2)

Mô hình converge sau vài phút và đây là kết quả phân bố hậu nghiệm và chuỗi MCMC của các hằng số trong mô hình. Ta quan tâm đến ka và ke. Giá trị của 2 hằng số này lần lượt là 1.5 và 0.08

plot(bpkmod, N = 5, ask = FALSE,theme=my_theme(8))

summary(bpkmod)

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: conc ~ Dose * ka * InV/(ka - ke) * (exp(-ke * Time) - exp(-ka * Time)) 
##          InV ~ 1 | Subject
##          ka ~ 1
##          ke ~ 1
##    Data: dat (Number of observations: 132) 
## Samples: 2 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 2000
## 
## Group-Level Effects: 
## ~Subject (Number of levels: 12) 
##                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(InV_Intercept)     0.49      0.15     0.28     0.86        178 1.00
## 
## Population-Level Effects: 
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## InV_Intercept     1.97      0.18     1.56     2.28        281 1.00
## ka_Intercept      1.52      0.13     1.28     1.80       1028 1.00
## ke_Intercept      0.08      0.01     0.07     0.09       1126 1.00
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma     1.16      0.08     1.02     1.31       1097 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Kết quả ka, ke này tương đương với cách làm cổ điển như trong website sau đây: http://sia.webpopix.org/nlme.html

Ta có thể kiểm tra xem mô hình có ước tính chính xác diễn tiến của biến conc theo thời gian hay không ? Trong hình dưới đây: Màu xanh chỉ giá trị thực tế, màu đỏ là giá trị mô hình ước tính.

pred_df=predict(bpkmod)%>%as.data.frame()

dat$Pred=pred_df$Estimate

post_df=bpkmod$fit%>%as.matrix()

ggplot(data=dat)+
  geom_point(aes(x=Time,y=conc),col="pink")+
  geom_line(aes(x=Time,y=conc),col="blue",linetype=2)+
  geom_line(aes(x=Time,y=Pred),col="red")+
  facet_wrap(~Subject,ncol=3,scales="free")+my_theme()

4 Kết luận

Package brms hỗ trợ thực hiện mô hình phi tuyến tính, phân cấp với cấu trúc mô hình có thể mở rộng thoải mái. Trong bài này chúng ta chỉ mới sử dụng mô hình 1 ngăn, và chỉ quan tâm đến hằng số; tuy nhiên các mô hình nhiều ngăn cũng có thể được giải quyết theo cùng cách làm này. Ta cũng có thể đưa thêm nhiều tham số vào mô hình và chúng có liên hệ với dữ liệu: thí dụ sử dụng cân nặng của bệnh nhân trong các mô hình dược động học, hay các chỉ số nhân trắc về hệ hô hấp cho các mô hình khí động học.

Vấn đề khó khăn nhất khi sử dụng phương pháp Bayes là xác định prior phù hợp. Với mô hình phi tuyến tính thì prior có vai trò cực kì quan trọng, như các bạn đã thấy trong bài.

Chúc các bạn thực hành vui. Hẹn gặp lại trong một bài khác về Bayes.