Cập nhật brms

Bài 1: Bayes Logistic GAM

1 Giới thiệu brms

Thân chào các bạn, đây là bài đầu tiên trong series cập nhật về phiên bản mới nhất của package brms. Cho những bạn chưa biết: brms là một giao thức mô hình hồi quy theo trường phái Bayes, thông qua ngôn ngữ STAN.

Có nghĩa là brms có họ hàng với những package sử dụng ngôn ngữ STAN khác như rstanarm, bayesplot, rstan… nhưng brms không tương đương với những package hồi quy Bayes khác sử dụng Gibs sampler (JAGs)- như MCMCglmm, và dĩ nhiên brms đối nghịch với toàn bộ phần còn lại của những package hồi quy khác trong R như glm đến gam, lme4, gamlss.

Như vậy brms không sử dụng quy tắc REML (ước tính hợp lý - likelihood cực đại) mà sử dụng phương pháp Bayes, nó ước tính phân phối hậu định của các tham số hồi quy trong mô hình, khi biết data, prior và likelihood function của model.

Ngay từ phiên bản đầu tiên ra mắt cách đây 2 năm, brms đã được cộng đồng tán thưởng vì tính giản dị của nó. Từ khi có brms, người dùng có thể thực hành hồi quy Bayes thoải mái mà không cần phải code mô hình trực tiếp bằng ngôn ngữ STAN.

Theo thời gian, brms ngày càng chứng tỏ ưu thế vượt trội so với rstanarm. Khả năng của brms ngày càng mở rộng. Vào tháng 8 năm 2017, brms đã được cải tiến vượt bậc với cấu trúc hoàn toàn mới để trở thành package hồi quy Bayes mạnh nhất, tương đương với gamlss bên phái REML.

Hiện nay brms giống như một con bạch tuộc mạnh mẽ và mềm dẻo:

brms có thể làm được:

1. Mô hình hóa cho từng tham số riêng biệt của các quy luật phân phối phổ biến (như trong GAMLSS), gồm Gaussian, student, skew_normal, Beta (liên tục), binomial, Bernoulli (nhị phân), Poisson, nhị thức âm, hurdle và các phân bố Zero/One inflated (số đếm), gamma.

2. Mô hình đa cấp, với random / grouping effects (như lme4 và lmer).

3. Sử dụng các hàm Smoothing để tạo ra mô hình generalized additive (GAM), bao gồm splines, tensor product, Gaussian process (như gam, mgcv).

4. Các mô hình phi tuyến tính phức tạp với nhiều parameter, thí dụ exponential.

5. Các hàm hiệu ứng đặc biệt, như monotonic, measurement error (bao gồm Meta-analysis).

6. Các mô hình Survival (Weibull, inverse.gaussian, lognormal) và Time series.

Để làm được những điều trên, cấu trúc của brms cũng như cú pháp các hàm đã thay đổi rất nhiều so với năm 2016. Ngoài ra một số package khác cũng được phát triển để hỗ trợ khai thác kết quả brms, bao gồm broom,bayesplot, sjPlot…

Vì những lý do này, nhóm BAV sẽ thực hiện một series bài thực hành hoàn toàn mới cho brms, để cập nhật những tính năng mới nhất của nó.

2 Mục tiêu

Trong bài đầu tiên này, chúng ta sẽ làm quen với một trong những tính năng mới của brms, đó là tích hợp những hàm Smoothing (splines, tensor product…) vào mô hình GLM và tạo ra mô hình GAM. Trong thí dụ minh họa sau đây, Nhi sẽ thực hiện một mô hình Bayesian Logistic GAM với phân bố Bernoulli có kèm tensor product function cho mỗi predictor.

Dataset Kyphosis (John M. Chambers, 1992) được chọn để minh họa cho thí nghiệm, đây là một data điển hình cho trường hợp quan hệ phi tuyến tính giữa Outcome và predictor. Dataset gồm 81 trẻ em với outcome là biến chứng Gù (Có/không) sau khi phẫu thuật cột sống. Hai predictor chính là Tuổi (tháng) và Vị trí đốt sống bị can thiệp (Start). Chúng ta sẽ giải quyết nó như một bài toán phân loại (mục tiêu xây dựng một quy tắc tiên lượng) và diễn dịch (giải thích quan hệ giữa predictor và nguy cơ xảy ra biến chứng gù hậu phẫu).

library(pander)
library(tidyverse)

df=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/rpart/kyphosis.csv")%>%as_tibble()

df=df%>%mutate(response=if_else(Kyphosis=="absent",0,1),
               Age=as.numeric(.$Age),
               Start=as.numeric(.$Start),
               Number=as.numeric(.$Number),
               X=as.factor(.$X))

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

X	Kyphosis	Age	Number	Start	response
1	absent	71	3	5	0
2	absent	158	3	14	0
3	present	128	4	5	1
4	absent	2	5	1	0
5	absent	1	4	15	0
6	absent	1	2	16	0

3 Thăm dò dữ liệu

Trước hết ta thăm dò dữ liệu bằng biểu đồ:

Boxplots cho thấy có sự tương phản khá rõ cho biến Start giữa 2 nhóm Có và không có Kyphosis, tuy nhiên phân phối của Start không đồng dạng, không đối xứng. Sự tương phản cho biến Age khó nhận thấy hơn.

df%>%gather(Age,Start,key="Features",value="Value")%>%
  ggplot(aes(x=Kyphosis,y=Value,fill=Kyphosis))+
  geom_boxplot(alpha=0.8)+
  theme_bw()+
  facet_wrap(~Features,scales="free",ncol=1)+coord_flip()

Khi xét Outcome như một biến số nhị phân với 2 giá trị 0/1 và dùng một mô hình binomial logistic đơn biến làm trung gian, ta có thể hình dung về quan hệ riêng phần giữa Age, Start và xác suất có Kyphosis:

binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), ...)
}

df%>%gather(Age,Start,key="Features",value="Value")%>%
  ggplot(aes(x=Value,y=response,col=Features,fill=Features))+
  geom_point()+
  binomial_smooth(formula = y ~ x,alpha=0.3)+
  theme_bw()+
  facet_wrap(~Features,scales="free",ncol=2)

4 Mô hình logistic GLM

Ta bắt đầu bằng một mô hình logistic cổ điển bằng hàm glm, phân phối binomial:

Kết quả mô hình tuyến tính cho thấy vị trí phẫu thuật (Start) là yếu tố duy nhất có ý nghĩa làm thay đổi nguy cơ bị kyphosis.

# logistic glm

lmod=glm(response~Age+Start,data=df,family=binomial(link="logit"))

summary(lmod)

## 
## Call:
## glm(formula = response ~ Age + Start, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5624  -0.5976  -0.3923  -0.2354   2.1998  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.224996   0.750638   0.300 0.764376    
## Age          0.009509   0.005904   1.611 0.107284    
## Start       -0.237939   0.064535  -3.687 0.000227 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 83.234  on 80  degrees of freedom
## Residual deviance: 65.299  on 78  degrees of freedom
## AIC: 71.299
## 
## Number of Fisher Scoring iterations: 5

pred0=predict(lmod,type="response")

preddf0=data_frame(Truth=df$Kyphosis,
                   Prob=pred0,
                   Predicted=if_else(pred0>0.5,"present","absent"))

caret::confusionMatrix(preddf0$Truth,preddf0$Predicted,positive="present")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction absent present
##    absent      59       5
##    present     11       6
##                                           
##                Accuracy : 0.8025          
##                  95% CI : (0.6991, 0.8827)
##     No Information Rate : 0.8642          
##     P-Value [Acc > NIR] : 0.9569          
##                                           
##                   Kappa : 0.3157          
##  Mcnemar's Test P-Value : 0.2113          
##                                           
##             Sensitivity : 0.54545         
##             Specificity : 0.84286         
##          Pos Pred Value : 0.35294         
##          Neg Pred Value : 0.92187         
##              Prevalence : 0.13580         
##          Detection Rate : 0.07407         
##    Detection Prevalence : 0.20988         
##       Balanced Accuracy : 0.69416         
##                                           
##        'Positive' Class : present         
## 

library(sjPlot)

sjt.glm(lmod,show.aic=T,
        show.hoslem = T,
        show.chi2 = T,
        show.r2=T,
        p.zero=T,
        emph.p=T,
        digits.p=6,
        digits.est = 3,
        digits.ci = 3)

		response
		Odds Ratio	CI	p
(Intercept)		1.252	0.281 – 5.614	0.764376
Age		1.010	0.998 – 1.022	0.107284
Start		0.788	0.687 – 0.888	0.000227
Observations		81
Pseudo-R2		R2CS = 0.199 R2N = 0.309 D = 0.221
AIC		71.299
Χ2deviance		p=0.000
Hosmer-Lemeshow-Χ2		5.127; p=0.744

sjp.glm(lmod)

sjp.glm(lmod,type="eff")

sjp.glm(lmod,type="ma")

## NULL

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: response
## 
## Terms added sequentially (first to last)
## 
## 
##       Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     80     83.234              
## Age    1    1.302        79     81.932    0.2539    
## Start  1   16.633        78     65.299 4.534e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sjp.glm(lmod,type="slope")

Tuy nhiên phẩm chất của mô hình logistic tuyến tính này còn rất hạn chế. Trên 2 biểu đồ thể hiện residuals theo 2 biến Age và Start, ta có thể dự đoán về một liên hệ phi tuyến tính giữa Outcome và predictors.

Để đưa hiệu ứng phi tuyến tính này vào mô hình, ta có thể dùng một số hàm trung gian, thí dụ: đa thức, smoothing spline, tensor product, Gaussian process… hoặc dùng một mô hình cây (decision tree với algorithm CART).

5 Decision tree

Ta thử mô hình cây trước:

library(rpart)
library(visNetwork)

tree<- rpart(Kyphosis ~ Age + Start, data = df)

tree

## n= 81 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 81 17 absent (0.79012346 0.20987654)  
##    2) Start>=8.5 62  6 absent (0.90322581 0.09677419)  
##      4) Start>=14.5 29  0 absent (1.00000000 0.00000000) *
##      5) Start< 14.5 33  6 absent (0.81818182 0.18181818)  
##       10) Age< 55 12  0 absent (1.00000000 0.00000000) *
##       11) Age>=55 21  6 absent (0.71428571 0.28571429)  
##         22) Age>=111 14  2 absent (0.85714286 0.14285714) *
##         23) Age< 111 7  3 present (0.42857143 0.57142857) *
##    3) Start< 8.5 19  8 present (0.42105263 0.57894737) *

plotcp(tree)

visNetwork::visTree(tree)

library(partykit)
plot(as.party(tree))

pred=predict(tree,df)

predtree=data_frame(Truth=df$Kyphosis,
                   Prob=pred[,2],
                   Predicted=if_else(Prob>0.5,"present","absent"))

caret::confusionMatrix(predtree$Truth,predtree$Predicted,positive="present")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction absent present
##    absent      53      11
##    present      2      15
##                                           
##                Accuracy : 0.8395          
##                  95% CI : (0.7412, 0.9117)
##     No Information Rate : 0.679           
##     P-Value [Acc > NIR] : 0.0008572       
##                                           
##                   Kappa : 0.5948          
##  Mcnemar's Test P-Value : 0.0265003       
##                                           
##             Sensitivity : 0.5769          
##             Specificity : 0.9636          
##          Pos Pred Value : 0.8824          
##          Neg Pred Value : 0.8281          
##              Prevalence : 0.3210          
##          Detection Rate : 0.1852          
##    Detection Prevalence : 0.2099          
##       Balanced Accuracy : 0.7703          
##                                           
##        'Positive' Class : present         
## 

Mô hình CART có thể đáp ứng cả 2 mục tiêu : tiên lượng và diễn dịch. Nó có độ chính xác khá cao so với logistic model (BAC=0.77)

Tuy nhiên vẫn còn một giải pháp khác, đó là GAM (generalized addtive model)

6 Mô hình GAM

Ta có thể áp dụng 2 hàm spline bậc 2 cho Age và cho Start, rồi thăm dò lại lần nữa trên biểu đồ logistic glm: Kết quả này có hình dạng khá giống với biểu đồ Residual vs predictor

df%>%gather(Age,Start,key="Features",value="Value")%>%
  ggplot(aes(x=Value,y=response,col=Features,fill=Features))+
  geom_point()+
  binomial_smooth(formula = y ~ splines::ns(x,2),alpha=0.3)+
  theme_bw()+
  facet_wrap(~Features,scales="free",ncol=2)

Bây giờ ta thử dựng một mô hình logistic GAM với package gam, sử dụng hàm spline bậc 2 cho Age và bậc 3 cho Start

Ghi chú: nhược điểm của mô hình GAM đó là predictor bị hoán chuyển thành kết quả của hàm Splines, nên rất khó diễn giải. ANOVA chỉ cho thấy ns(Age,2) và s(Start,3) gây hiệu ứng ý nghĩa làm thay đổi Outcome.

library(gam)

fitgam=gam(response ~ ns(Age,2) + s(Start,3), data=df, 
           family=binomial)

summary(fitgam)

## 
## Call: gam(formula = response ~ ns(Age, 2) + s(Start, 3), family = binomial, 
##     data = df)
## Deviance Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.67213 -0.43435 -0.24601 -0.04133  2.15865 
## 
## (Dispersion Parameter for binomial family taken to be 1)
## 
##     Null Deviance: 83.2345 on 80 degrees of freedom
## Residual Deviance: 50.5232 on 75 degrees of freedom
## AIC: 62.5233 
## 
## Number of Local Scoring Iterations: 9 
## 
## Anova for Parametric Effects
##             Df Sum Sq Mean Sq F value   Pr(>F)   
## ns(Age, 2)   2  5.906  2.9529  4.3865 0.015791 * 
## s(Start, 3)  1  6.577  6.5774  9.7705 0.002523 **
## Residuals   75 50.490  0.6732                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar Chisq  P(Chi)  
## (Intercept)                             
## ns(Age, 2)                              
## s(Start, 3)       2     6.9065 0.03164 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Một cách làm khác cho GAM là dùng package gmcv:

Đây là mặt phẳng hồi quy của một mô hình logistic tuyến tính (như glm)

library(mgcv)

fit <- gam(response~Age+Start,
            family=binomial(),
            data=df)

vis.gam(fit,view=c("Age","Start"),theta=50,phi=50,n.grid=30,color="topo",ticktype = "detailed") 

Khi áp dụng 2 hàm splines cho Age và Start, mặt phẳng hồi quy trong không gian thay đổi như sau:

fit2 <- gam(response~s(Age)+s(Start),
           family=binomial(),
           data=df)

vis.gam(fit2,view=c("Age","Start"),theta=35,phi=30,n.grid=50,color="topo",ticktype = "detailed") 

summary(fit2)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## response ~ s(Age) + s(Start)
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.1723     0.4846  -4.483 7.37e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##            edf Ref.df Chi.sq p-value   
## s(Age)   2.183  2.761  6.346 0.08054 . 
## s(Start) 2.199  2.740 12.798 0.00421 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =   0.33   Deviance explained = 37.1%
## UBRE = -0.2204  Scale est. = 1         n = 81

7 Bayes logistic GAM trong brms

Bây giờ chúng ta đi vào phần chính, đó là mô hình Bayesian GAM sử dụng brms.

Một số chi tiết cần giải thích như sau:

1. brms sử dụng một giao thức đơn giản, cho phép người dùng mô tả nội dung model bằng formula như trong glm hay lme4.

2. Nội dung này sau đó sẽ được phiên dịch thành STAN code, với một prior mặc định (người dùng có thể tùy chỉnh prior).

3. Brms xử lý data dưới dạng vector cho Outcome và Matrix cho predictors.

4. Khác với glm, brms áp dụng phân phối Bernoulli với link function là logit cho logistic model chứ không phải là phân phối Binomial (dĩ nhiên bạn cũng có thể dùng Binomial, nhưng STAN sẽ chạy chậm hơn).

5. Trong STAN code, brms cũng tính cả log_likelihood cho model, từ đó cho phép tính WAIC và LOOIC.

6. Phần thứ 2 trong syntax brms là điều kiện lấy mẫu cho những chuỗi MCMC, với quy tắc hoàn toàn giống như trong rstan.

7. Khi thi hành mô hình, đầu tiên STAN code sẽ được phiên dịch thành một chương trình sử dụng ngôn ngữ C++, và nếu không có lỗi syntax, sampler sẽ bắt đầu chạy để tạo các chuỗi MCMC.

8. Sau khi converge, ta sẽ có kết quả là một mô hình STAN chứa phân phối hậu nghiệm của các tham số.

Cách diễn giải và khai thác mô hình này đã được Nhi trình bày nhiều lần trong các tutorial về Bayes.

Cấu trúc mô hình Bayes trong brms như sau:

Mục tiêu của mô hình là ước tính một biến kết quả Y từ hàm mật độ xác suất của một quy luật phân phối D, được xác định từ những tham số theta:

\[y \sim D(\theta _{1},\theta _{2},...)\]

Mỗi tham số theta (thí dụ Mu, Sigma, Nu, Zero_inflated,…) có thể được xem như một outcome thứ cấp, và được mô hình hóa bởi một link function fp cho một trị số eta(p), thí dụ trong mô hình logistic thì fp=logit và phân phối D = binomial

\[\theta _{p} \sim f_{p}(\eta_{pi})\]

Eta được xem là outcome của một mô hình tuyến tính có công thức:

\[\eta = \beta X + Zu + \sum_{K}^{1}s_{k}(x_{k})\]

Trong đó X là model matrix cho population hay fixed effect (hiệu ứng chính), beta là các tham số hồi quy cho matrix X này.

Mô hình này có thể tích hợp thêm: Z là Hiệu ứng ngẫu nhiên do phân nhóm (random effects, grouping effects), u là tham số của Z; một hay nhiều additive terms bậc K cho một predictor, bao gồm hàm Smoothing splines (s), Gaussian process (gp) hay Tensor product (t2).

brms mượn phần lõi của nlme để ước tính random effect và các hàm splines trong mgcv để tạo mô hình GAM.

Nhi sử dụng hàm t2 (tensor product) thay cho hàm spline, bạn cũng có thể dùng s() nữa, tốc độ converge là như nhau. hàm gp() (Gaussian process) chạy chậm hơn nhiều và kết quả không khả quan bằng.

Lưu ý là model này chưa có prior

set.seed(123)

library(brms)

bmod <- brm(response ~ t2(Age)+t2(Start), 
            data=df,family=bernoulli(link="logit"),
                 chains = 2,
            iter=2000,
            warmup=500,
            cores = 4,
            control=list(adapt_delta=0.99))

8 Khai thác mô hình Bayes

Kết quả của mô hình STAN có thể được khai thác khá phong phú như sau:

1. Đầu tiên ta có thể trích xuất STAN code cho mô hình:

bmod$model

## // generated with brms 1.10.2
## functions { 
## } 
## data { 
##   int<lower=1> N;  // total number of observations 
##   int Y[N];  // response variable 
##   int<lower=1> K;  // number of population-level effects 
##   matrix[N, K] X;  // population-level design matrix 
##   // data of smooth t2(Age)
##   int nb_1;  // number of bases 
##   int knots_1[nb_1]; 
##   matrix[N, knots_1[1]] Zs_1_1; 
##   // data of smooth t2(Start)
##   int nb_2;  // number of bases 
##   int knots_2[nb_2]; 
##   matrix[N, knots_2[1]] Zs_2_1; 
##   int prior_only;  // should the likelihood be ignored? 
## } 
## transformed data { 
##   int Kc = K - 1; 
##   matrix[N, K - 1] Xc;  // centered version of X 
##   vector[K - 1] means_X;  // column means of X before centering 
##   for (i in 2:K) { 
##     means_X[i - 1] = mean(X[, i]); 
##     Xc[, i - 1] = X[, i] - means_X[i - 1]; 
##   } 
## } 
## parameters { 
##   vector[Kc] b;  // population-level effects 
##   real temp_Intercept;  // temporary intercept 
##   // parameters of smooth t2(Age)
##   vector[knots_1[1]] zs_1_1; 
##   real<lower=0> sds_1_1; 
##   // parameters of smooth t2(Start)
##   vector[knots_2[1]] zs_2_1; 
##   real<lower=0> sds_2_1; 
## } 
## transformed parameters { 
##   vector[knots_1[1]] s_1_1 = sds_1_1 * zs_1_1; 
##   vector[knots_2[1]] s_2_1 = sds_2_1 * zs_2_1; 
## } 
## model { 
##   vector[N] mu = Xc * b + Zs_1_1 * s_1_1 + Zs_2_1 * s_2_1 + temp_Intercept; 
##   // priors including all constants 
##   target += normal_lpdf(zs_1_1 | 0, 1); 
##   target += student_t_lpdf(sds_1_1 | 3, 0, 10)
##     - 1 * student_t_lccdf(0 | 3, 0, 10); 
##   target += normal_lpdf(zs_2_1 | 0, 1); 
##   target += student_t_lpdf(sds_2_1 | 3, 0, 10)
##     - 1 * student_t_lccdf(0 | 3, 0, 10); 
##   // likelihood including all constants 
##   if (!prior_only) { 
##     target += bernoulli_logit_lpmf(Y | mu); 
##   } 
## } 
## generated quantities { 
##   // actual population-level intercept 
##   real b_Intercept = temp_Intercept - dot_product(means_X, b); 
## }

2. Hàm summary để xem nội dung mô hình:

summary(bmod)

##  Family: bernoulli(logit) 
## Formula: response ~ t2(Age) + t2(Start) 
##    Data: df (Number of observations: 81) 
## Samples: 2 chains, each with iter = 2000; warmup = 500; thin = 1; 
##          total post-warmup samples = 3000
##     ICs: LOO = NA; WAIC = NA; R2 = NA
##  
## Smooth Terms: 
##                Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sds(t2Age_1)       5.84      3.92     0.87    16.18       1404 1.00
## sds(t2Start_1)     6.10      4.93     0.80    19.10       1130 1.00
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept    -3.85      1.24    -6.93    -1.99        913 1.00
## t2Age_1      -0.06      0.97    -1.67     2.25       1350 1.00
## t2Start_1    -1.98      0.85    -4.11    -0.76        880 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).

Nội dung model Bayes GAM phân biệt 2 nhóm : Smooth terms (sự tham gia của hàm tensor product) và Population-Level effects (hiệu ứng chính).

Smooth terms cho phép kết luận về quan hệ phi tuyến tính giữa outcome Eta và predictor thông qua hàm smoothing.

Trong khi đó phần hiệu ứng chính còn lại cho thấy:

Biến t2_Age không có liên hệ ý nghĩa với kyphosis (95%CI chứa giá trị 0), còn giá trị t2_ Start nghịch biến với nguy cơ kyphosis, tức là Start thấp hơn tương ứng với nguy cơ cao hơn.

Kết quả phân phối hậu định của smoothing term và fixed effect có thể được tóm tắt bằng hàm tidy của package broom:

library(broom)

bmod %>% tidy(parameters = c("^sds_","^b_"),prob = 0.975)

##            term    estimate std.error      lower      upper
## 1   b_Intercept -3.84984033 1.2431667 -7.5683631 -1.8302751
## 2     b_t2Age_1 -0.05917458 0.9716609 -1.8439455  2.8952187
## 3   b_t2Start_1 -1.98410535 0.8462535 -4.7672140 -0.6502784
## 4   sds_t2Age_1  5.84452998 3.9246924  0.4164998 19.1748280
## 5 sds_t2Start_1  6.09753661 4.9251792  0.3746217 23.2810040

3. Ta có thể kiểm tra đặc tính mô hình Bayes bằng tiêu chuẩn LOOIC, WAIC và R2

loo(bmod)

##  LOOIC    SE
##  67.41 11.83

waic(bmod)

##   WAIC    SE
##  66.53 11.56

bayes_R2(bmod)

##    Estimate  Est.Error  2.5%ile  97.5%ile
## R2 0.358799 0.07382833 0.193539 0.4766589

4. Về biểu đồ, ta có nhiều dạng như: marginal effect plot

marginal_effects(bmod)

5. Hàm plot sẽ vẽ phân phối hậu định và các chuỗi MCMC

bmod%>%plot(combo = c("hist", "trace"), widths = c(1,1.5),
  theme = theme_bw(7))

6. Đặc biệt package sjPlot cũng áp dụng được cho mô hình STAN, và lần này nó vẽ ra phân phối liên tục của Odds-ratio luôn cho chúng ta :

plot_model(bmod,transform = "exp",bpe.style = "dot",colors="Set1")

7. Ta có thể kiểm tra một giả thuyết H1 tùy thích cho Odds-ratio bằng Bayes Factor, sử dụng hàm hypothesis

Thí dụ nếu ta muốn kiểm tra H1= Odds-ratio của t2(Start) nhỏ hơn 0.5, kết quả là Bayes Factor cao hơn 57.

Nếu ta tăng ngưỡng H1 này lên 0.8, Bayes Factor sẽ cao gần 1500 :

hypothesis(bmod,"exp(t2Start_1)<0.5")

## Hypothesis Tests for class b:
##                          Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(t2Start_1))-... < 0    -0.32      0.12     -Inf     -0.1      57.82
##                          Star
## (exp(t2Start_1))-... < 0    *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

hypothesis(bmod,"exp(t2Start_1)<0.8")

## Hypothesis Tests for class b:
##                          Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(t2Start_1))-... < 0    -0.62      0.12     -Inf     -0.4       1499
##                          Star
## (exp(t2Start_1))-... < 0    *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

Ngược lại BF cho giả thuyết OR của Age <1 là rất thấp

hypothesis(bmod,"exp(t2Age_1)<1")

## Hypothesis Tests for class b:
##                        Estimate Est.Error l-95% CI u-95% CI Evid.Ratio
## (exp(t2Age_1))-(1) < 0     0.93      8.42     -Inf     4.49       1.39
##                        Star
## (exp(t2Age_1))-(1) < 0     
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

8. Cuối cùng, ta có thể sử dụng mô hình Bayes logistic GAM này để tiên lượng :

Mô hình này có phẩm chất tốt hơn nhiều so với mô hình logistic tuyến tính : nó có AUC = 89 và BAC=0.79

pred=predict(bmod,df)%>%.[,1]

library(classifierplots)

roc_plot(df$response, pred)+theme_bw()

## [1] "Calculating AUC ..."
## [1] "(AUC) Sorting data ..."
## [1] "(AUC) Calculating ranks ..."
## [1] "AUC: 89.3382352941177"
## [1] "Bootstrapping ROC curves"
## [1] "Eval AUC"
## [1] "Producing ROC plot"

preddf=data_frame(Truth=df$Kyphosis,
                  Prob=pred,
                  Predicted=if_else(pred>0.5,"present","absent"))

caret::confusionMatrix(preddf$Truth,preddf$Predicted,positive="present")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction absent present
##    absent      60       4
##    present      8       9
##                                          
##                Accuracy : 0.8519         
##                  95% CI : (0.7555, 0.921)
##     No Information Rate : 0.8395         
##     P-Value [Acc > NIR] : 0.4534         
##                                          
##                   Kappa : 0.5111         
##  Mcnemar's Test P-Value : 0.3865         
##                                          
##             Sensitivity : 0.6923         
##             Specificity : 0.8824         
##          Pos Pred Value : 0.5294         
##          Neg Pred Value : 0.9375         
##              Prevalence : 0.1605         
##          Detection Rate : 0.1111         
##    Detection Prevalence : 0.2099         
##       Balanced Accuracy : 0.7873         
##                                          
##        'Positive' Class : present        
## 

9 Tổng kết

Nhi vừa trình bày một trong số những tính năng mới mà brms hỗ trợ, đó là mô hình GAM với smoothing spline. Cách làm này đã cải thiện đáng kể tính chính xác của mô hình cho mục tiêu tiên lượng.

Khi sử dụng brms, các bạn đang thực hành hồi quy Bayes một cách thoải mái mà không hề hay biết. Khác biệt duy nhất nằm ở phần diễn giải và khai thác kết quả.

Trong những bài tiếp theo, chúng ta sẽ trải nghiệm những tính năng mạnh mẽ và linh hoạt khác mà brms có thể làm được. Hy vọng các bạn có thể tự tin hơn khi chọn đi theo trường phái Bayes.